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1.     INTRODUCTION 

The  purpose  of  this  report  is  to  discover  and  quantify  the  systematic 
errors  in  the  algorithms  employed  by  NUWES  in  their  short  base  line 
underwater  position  location  system.  Systematic  error  can  have  a  number  of 
sources  and  previous  works  [6,7]  have  treated  other  issues.  The  present  work 
deals  with  the  algorithms  used  in  sound  ray  tracing.  There  are  a  number  of 
aspects  to  be  treated,  and  some  background  is  necessary  in  order  to  explain 
them. 

Figure  1  contains  a  schematic  diagram  of  a  short  base  line  hydrophonic 
array  and  of  the  signals  it  may  receive.  Sharply  pulsed  signals,  or  pings,  are 
sent  by  the  sound  source  vehicles  (surface  craft,  submarine,  torpedo)  and  they 
are  received  by  the  four  transducers  (called  the  X,  Y,  Z  and  C-phones)  of  the 
array.  These  four  hydrophones  form  a  right  angled  coordinate  system  with 
origins  at  the  C-phone  and  arm  lengths  D  (=30  feet)  to  each  of  the  other  three 
phones. 

The  ray  paths  from  a  source  to  the  four  receivers  are  synchronously  timed 
with  great  precision  (10  sees)  and  the  differentials  of  arrival  times  are  used 
to  construct  the  direction  of  the  source.  But  due  to  variability  of  the  speed  of 
sound  at  various  water  depths,  the  ray  paths  themselves  are  not  straight  lines. 
Also  the  paths  may  change  from  day  to  day  as  the  water  depth-velocity  profile 
changes.  Knowledge  of  the  current  profile  allows  one  to  perform  a  ray  tracing 
computation.  Recovery  of  the  three-dimensional  position  of  the  sound 
source  is  accomplished  by  reconstructing  the  ray  path  and  following  it  for  the 
given  amount  of  time.  Each  source  vehicle  has  a  phase  coded  "ping"  so  that 
its  signals  can  be  discriminated  from  those  of  other  vehicles. 
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Figure  1.  Short  Baseline  Array  and  Signal  Sources 


Each  hydrophonic  array  is  placed  on  the  sea  floor  by  lowering  it  over  the 
side  of  a  utility  ship.  Each  has  a  self-leveling  capability.  When  finally  resting 
on  the  bottom,  they  are  not  perfectly  level,  but  the  X  and  Y  arms  have  tilt 
meters  that  measure  the  angles  that  are  made  with  the  horizontal.  An  array 
surveying  activity  checks  these  angles  and  measures  the  rotation  of  the 
vertical.  Thus  the  local  coordinate  system  can  be  reconciled  with  the  master 
range  coordinate  system. 

The  currently  employed  ray  tracing  methodology  partitions  the  water 
column  into  a  number  of  equal  thickness  layers  and  treats  the  speed  of  sound 
as  constant  in  each  layer.  This  leads  to  the  use  of  isospeed  ray  tracing  [1,  5]. 
Presently  the  layer  thicknesses  at  Nanoose  are  25  feet.  In  those  instances  for 
which  the  measured  depth-velocity  profile  does  not  go  as  deep  as  the  array,  a 
constant  speed  extrapolation  is  used.  In  effect  the  thickness  of  the  deepest 
layer  is  larger,  perhaps  50  or  75  feet. 

Now  the  issues  can  be  detailed: 

i)       How  accurate  is  iso-speed  ray  tracing  using  a  25  foot  water  layer 
increment? 

ii)     What  is  the  effect  when  it  is  necessary  to  use  a  thicker  deepest 
layer? 

iii)    How  well  are  the  ray  tracing  directions  and   transit   times 
determined? 

iv)    What  effect  does  the  various  depth  velocity  profiles  have  upon 
the  answers  to  i),  ii)  and  iii)? 

The  treatment  of  these  questions  requires  a  valid  sound  ray  construction 

methodology  and  a  representative  set  of  depth-velocity  profiles.   To  satisfy  the 

latter  requirement  we  have  selected  twelve  experimental  days  at  the  Nanoose 

range  spanning  the  period  May  1988  to  June  1989.    They  are  presented  (in 

graphical  form)  in  Appendix  A.    We  note  that  the  more  interesting  ones  seem 


to  appear  in  the  Spring.  Two  of  the  experimental  days,  23-24  April,  1989  are 
consecutive.  This  allows  us  a  glimpse  into  the  question  of  day  to  day 
variability. 

To  treat  the  former  requirement,  we  have  developed  an  isogradient  ray 
fitting  algorithm.  It  fits  a  sound  ray  connecting  two  given  points  (in  the 
horizontal-vertical  plane)  assuming  direct  path  propagation.  The  depth 
velocity  profile  is  partitioned  into  five  foot  equi  thickness  layers  and  within 
each  layer  the  slope  of  sound  speed  vs.  depth  is  constant.  Thus  the  DV  profile 
is  represented  as  a  continuous  function  consisting  of  a  sequence  of  straight 
line  segments.  Within  each  layer  the  ray  path  is  a  circle  arc  because  of  Snell's 
law,  [2,3,8].  The  outputs  of  this  algorithm  are  the  angles  that  the  ray  makes 
with  the  horizontal  at  each  of  the  endpoints,  and  the  transit  time  of  the  ray 
from  the  initial  point  to  the  final  one. 

Section  2  of  the  report  presents  some  theoretical  material  and  formulas. 
The  distinctions  between  isospeed  and  isogradient  ray  tracing  are  explained. 
Snell's  law  is  introduced  and  supported. 

Section  3  of  the  report  deals  with  the  accuracy  of  pure  ray  tracing  in  two 
dimensions,  the  horizontal-vertical  plane.  Computations  are  made  for  a 
number  of  initial  angle-transit  time  pairs  and  for  all  twelve  depth  velocity 
profiles.  Errors  from  this  source  are  generally  small  but  can  be  as  large  as  a 
foot  or  more.  The  situation  is  more  difficult  when  extrapolation  of  the  water 
column  is  necessary.  This  occurs  when  the  sounding  does  not  extend  as  deep 
as  the  receiver.  In  such  cases  we  cannot  be  definite  about  the  nature  of  the 
errors,  but  several  equally  defensible  methods  lead  to  results  that  disagree  by 
five  or  ten  feet  and  even  more.     We  state  that  there  is  a  problem  with 


extrapolation.  The  soundings  should  be  made  at  the  deepest  part  of  the  range 
and  to  full  depth.  Failing  that,  a  careful  development  of  an  extrapolation 
policy  should  be  made. 

Section  4  deals  with  the  three  dimensional  problem  of  locating  the 
position  of  the  source  based  upon  the  transit  times  to  each  of  the  four 
hydrophones.  Hence  the  question  is  one  of  finding  the  azimuth  and  initial 
elevation  angles  of  the  ray,  and  a  matching  transit  time  to  stop  the  ray  tracing 
algorithm.  Also  there  are  confounding  sources  of  variability.  The  depth- 
velocity  profile  plays  a  role,  as  mentioned  before.  But  the  arrays  themselves 
have  directional  properties  that  would  interact  with  the  algorithm  even  if 
they  were  fully  level  and  aligned  with  the  range.  The  fact  that  the  arrays  are 
tilted  and  rotated  in  a  variety  of  ways  has  contributed  to  the  puzzle  of 
interpreting  the  mismatches  in  the  array  overlap  areas.  Taken  altogether  it  is 
shown  that  systematic  error  from  these  sources  can  be  as  large  as  ten  or 
twelve  feet.  Moreover  errors  of  this  magnitude  are  unnecessary.  An 
alternative  method  is  proposed  which  can  reduce  them  by  at  least  an  order  of 
magnitude. 

The  conclusions  are  summarized  in  Section  5.  Section  6,  an  addendum, 
addresses  an  issue  raised  in  the  review  process.  Also  a  number  of  appendices 
are  included.  They  hold  the  depth-velocity  profiles,  supporting  mathematical 
details,  details  of  the  algorithms,  and  the  source  code  for  the  FORTRAN 
programs. 

A  brief  general  statement  of  conclusions  is  as  follows: 

i)       The  error  in  iso-speed  raytracing  is  an  increasing  function  of 
horizontal  range,  but  is  seldom  more  than  one  foot. 


ii)  The  error  due  to  constant  speed  extrapolation  in  the  deepest 
layer  can  range  up  to  10  or  more  feet. 

iii)  The  error  due  to  initializing  the  ray  tracing  is  a  periodic 
function  of  the  azimuth  direction  and  can  be  substantial  for 
tilted  arrays  and  at  the  greater  horizontal  ranges.  The  effect  of 
the  determination  of  azimuth  is  especially  noticeable. 

Greater  details  are  presented  as  the  various  issues  are  developed. 

2.     PERTINENT  ITEMS  FROM  RAY  TRACING 

The  sound  source  pings  and  sends  out  an  isotropic  wave  front,  which  is  a 
fixed  phase  point  on  the  pressure  cycle.  A  receiving  transducer  measures  the 
time  of  arrival  of  the  wave  front.  A  ray  is  a  path  normal  to  the  wave  front 
and  extending  in  space  back  to  the  source.  Our  goal  is  to  trace  a  ray  from  the 
receiver  for  a  fixed  amount  of  time  and  thereby  locate  the  source.  To  do  this 
we  must  construct  the  azimuth  and  elevation  angles  of  the  ray  at  the  receiver 
and  then  bend  the  ray  back  through  the  various  speed  layers  until  the 
measured  transit  time  is  consumed. 

If  the  speed  of  sound  in  water  were  constant  then  the  ray  path  would  be  a 
straight  line.  But  it  is  not.  Speed  is  a  function  of  temperature,  pressure  and 
salinity.  These  variables  interact  in  interesting  ways  for  the  water  layer  that 
affects  our  problem.  Speed  is  not  necessarily  a  monotone  function  of  depth. 
Conditions  change  with  time,  and  water  sounding  drops  are  made  daily. 
They  provide  a  depth-velocity  profile  which  is  assumed  to  be  fixed  for  the 
entire  day's  exercises.  Further,  these  values  are  assumed  constant  throughout 
each  horizontal  plane;  i.e.,  the  field  is  homogeneous. 

Our  immediate  goal  is  to  justify  the  use  of  a  ray  invariant  in  a 
horizontal-vertical  plane  and  to  establish  the  circular  arc  nature  of  ray  paths 


in  water  layers  for  which  the  speed  of  sound  is  a  straight  line  function  of 
depth. 

Let    us    begin     with         a+b 
Snell's  law.    Consider  two  ^<T\ 

adjacent  layers  with  speed 
S]  in  the  lower  layer  and 
speed  S2  in  the  upper.  The 

ray  enters  the  lower  layer  at  elevation  angle  0i  and  the  upper  layer  with  angle 
02-  Given  si  and  S2  let  us  find  the  relationship  between  61  and  02  that  will 
minimize  the  transit  time  from  (0,0)  to  (c,a+b). 

Proposition.    For  a  ray  to  traverse  from  (0,0)  to  (c,a+b)  in  minimum  time, 
we  must  have  the  relationship 

cos(0i)     cos(02) 


si 


S2 


(2.1) 


Proof.    The  transit  time  of  a  path  from  (0,0)  to  (c,a+b)  that  goes  through 
(x,a)  is  given  by 

1 


T(x)  =  -  Va2  +  x2  +  -  Vb2  +  (c-x)2. 


Further,  it  has  derivative 


T'(x)  = 


1 


c-x 


si  Va2  +  x2      s2  Vb2  +  (c-x)2 


The  relationship  (2.1)  is  a  consequent  of  setting  T'(x)  =  0. 

Now  suppose  the  point  c  is  not  fixed  but  variable.    It  follows  that  a  ray 
entering  the  lower  layer  at  an  angle  0]  (with  the  horizontal)  will  seek  the  path 


of  minimum  transit  time  and  exit  the  upper  layer  at  an  angle  of  62;  the  two 
angles  are  related  by  (2.1). 

Next  suppose  that  a  number  of  layers  are  stacked  vertically;  within  each 
layer  the  speed  of  sound  is  constant.  The  relationship  (2.1)  must  hold  for 
every  successive  pair  of  layers  and  hence  the  ratio 

cos(6) 
rv  =  — —  (2.2) 

must  be  constant  for  the  entire  ray  path.   This  value  characterizes  the  ray  path 
and  is  called  the  ray  invariant.   This  is  Snell's  law. 

Consider  the  vertical  plane  containing  the  source  and  the  receiver.  Call 
this  is  (h,z)  plane  with  the  depth  z  taken  as  positive  downward.  Our  ray 
tracing  problem  is  two  dimensional  in  this  plane.  Given  the  depth  velocity 
profile  we  need  only  the  elevation  angle  at  the  receiver  and  the  transit  time 
to  locate  the  source. 

The  depth- velocity  profile  can  be  approximated  with  a  series  of  thin  layers 
each  having  constant  (internal)  sound  speed.  Thus  ray  tracing  can  be  enacted 
using  such  an  approximation.  This  approach  is  called  isospeed  ray  tracing  [1]. 

A  sharper  approximation  is  available  through  the  use  of  water  layers 
whose  sound  speed  structure  can  be  represented  with  a  linear  function  of 
depth 

u(z)  =  uo  +  u1z  (2.3) 


8 


Proposition.  If  (2.3)  holds  then 
the  ray  path  is  a  circle  arc  of 
radius  zm  +  Vo/\)\.  If  U]  >  0  then 
the  arc  is  a  distance  \)o/"Wi  above 
z  =  0. 

Proof.    Referring  to  the  diagram, 
let  (h,z)  be  an  arbitrary  point  on 
the  ray  path  and  0  is  the  angle  of  the  tangent  line. 
Because  of  Snell's  law  we  must  have 

cos(6)  1 


DO+UiZ      X)Q+  t>izmax 


=  rv 


and  the  ray  path  can  be  described  parametrically  in  terms  of 

{h(0),  z(9)}       for  0  <  0 


The  radius  of  curvation  can  be  found  from  the  general  formula 

3/2 
{[h'(6)]2  +  [Z'(0)]2} 

r  = 


dh 


h'(0)z"(0)  -  h"(0)z'(0)  | 

Using  (2.3),  -j—  =  cot(0),  and  implicit  differentiation  we  find 

z'(0)  =  -sin(0)/u!rv  h'(0)  =  -cos^/^rv 


z"(0)  =  -cos(0)/\)irv 


and  hence 


h"(0)  =  sin(0)/\)!rv 


uo 


r  =  l/\)irv=— +  zmax 


(2.4) 


which  is  a  constant,  independent  of  0.  Hence  the  path  is  a  circle  arc;  the 
radius  is  as  specified;  and  the  center  of  the  circle  is  on  a  line  above  zero  at  a 
distance  \)q/\)-[.  This  proof  follows  those  found  in  [4,8]. 

If  \)i  <  0,  then  the  [h(0),  z(0)]  curve  is  concave,  still  a  circle  arc,  and  the 
circle's  center  is  still  on  the  line  z  =  -uo/ui.  But  now  this  line  is  below  z=0.  If 
\>i  =  0,  the  circle  radius  is  infinite,  the  sound  speed  is  constant,  and  the  ray 
path  is  a  straight  line. 

Now,  the  numerical  construction  of  the  ray  path  can  also  be  accomplished 
by  representing  each  layer's  sound  speed  as  a  straight  line  segment  (function 
of  depth)  and  piecing  together  the  consequent  circle  arcs.  Since  each  layer  has 
a  constant  speed  gradient  with  depth,  this  is  called  isogradient  ray  tracing,  [1]. 

The  question  of  efficiency  of  the  two  methods,  isospeed  and  isogradient 
ray  tracing,  is  really  a  question  of  how  well  the  depth-speed  function  in  a  layer 
can  be  represented  by  a  constant  on  the  one  hand,  or  by  a  linear  function  on 
the  other.  For  thick  layers  there  may  be  oscillations  that  make  the  choice 
difficult.  For  thin  layers  it  seems  that  the  straight  line  fit  should  perform 
better. 

The  algorithm  for  isogradient  ray  tracing  is  presented  in  Appendix  C.  A 
corresponding  algorithm  for  isospeed  ray  tracing  can  be  extracted  from  it  by 
making  a  number  of  deletions.  Fortran  source  codes  for  each  are  shown  in 
Appendix  G.  Inputs  for  these  algorithms  include  the  depth-velocity  table;  the 
layer  depths;  the  depth  of  the  receiver;  the  elevation  (layer  entrance)  angle  at 
the  receiver;  the  ray  transit  time.  The  outputs  are  the  horizontal  and  vertical 
end  points  of  the  ray,  and  the  final  elevation  (exit  layer)  angle  of  the  ray.  This 
latter  quantity  is  needed  in  the  timing  synchronization  model,  [7]. 
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For  isospeed  ray  tracing  the  pertinent  formulae  are 

Ah  =  Az  cot(9) 

At  =  Az/\)  sin(0) 
where  0  is  the  angle  that  the  ray  enters  the  layer;  Az  is  the  layer  thickness;  Ah 
is  the  horizontal  distance  transversed  in  the  layer;  At  is  the  transit  time 
through  the  layer.    The  next  layer  entrance  angle  is  computed  from  the  ray 
invariant  equation  (2.1). 

For  isogradient  tracing  the  pertinent  formulas  are  more  complicated.  For 
a  ray  that  enters  a  layer  at  depth  zo;  horizontal  displacement  ho;  and  angle  80 
we  must  first  compute  the  coordinate  (qi,  qo)  of  the  center  of  the  circle  arc 
0)1*0) 


qi  =  -i)o /vi 

qi  =  ho  +  (q2-zo)  sin(0o)/cos(0o) 
and  the  radius  of  the  arc 

r  =  signum  (q2)  (q2  -  zo)/cos(0o). 

The  new  horizontal  displacement  is 

hi  =  qi  -  signum(q2)  r  sin(0o) 

and  the  increase  in  transit  time  is  the  line  integral  At  =  J 


(2.5) 


(2.6) 


(2.7) 


ds 


v0  +  v^s) 


along  the 


circular  path.   This  integral  is  most  easily  managed  using  ds  =  y{dz)    +(dh)     = 
dO  I  \)\rv  and  uo  +  Diz  =  cos(0)/rv  so  that 


Ol 


«=-hl 


dd 


vi  J  cos(^)      vi 


In 


1  +  sinf^) 
cos(0i) 


In 


l  +  sin(#o) 
cos(e0) 
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and  the  layer  exit  angle  is  computed  from  the  ray  invariant  equation  2.1), 

cos(0i)  =  rv\)(zi) 

The  layer  exit  angle  is  the  entrance  angle  for  the  next  layer. 

The  above  equations  form  the  heart  of  direct  path  ray  tracing.  The 
organizational  questions  that  arise  when  developing  a  ray  tracing  algorithm 
are  treated  in  Appendix  C.  Basically  one  proceeds  upwards  through  the  layers 
until  the  specified  total  transit  time  is  consumed.  An  end  correction  is 
normally  necessary  because  of  the  requirement  to  stop  part  way  through  a 
layer. 

3.     QUALITY  OF  ISOSPEED  RAY  TRACING 

We  are  concerned  with  the  quality  of  the  currently  employed  isospeed  ray 
tracing  algorithm,  which  uses  a  uniform  water  layer  thickness  of  25  feet.  The 
receiver  depths  range  from  about  1100  to  1350  feet.  The  elevation  angle  can 
range  from  90°  (directly  overhead)  to  some  rather  small  but  positive  values. 
Of  course  a  variety  of  water  columns  (depth-velocity  profiles)  can  be 
encountered.  We  have  selected  twelve  (see  Appendix  A)  spanning  the  period 
May,  1988  to  July  1989. 

Our  first  problem  is  to  establish  a  standard  ray  to  serve  as  a  basis  for 

comparison.     To  this  end  we  are  limited  by  the  resolution  of  the  water 

column  data  available.    The  information  depicted  graphically  in  Appendix  A 

provides  sound  speed  averages  for  every  five  feet.    That  is,  at  level  /  the 

corresponding  velocity  value  represents 

1+5 
vi  =-  \v{zyiz  (3.1) 

i 
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and  we  have  no  information  concerning  the  amount  of  variability  that  may 
exist  within  the  layer.  It  is  presumed  small  and  is  neglected.  (A  model  for 
assessing  such  variability  is  presented  in  Appendix  D,  and  this  author  is 
concerned  about  the  issue  for  small  entrance  angles.) 

The  most  expedient  standard  available  is  to  employ  the  ray  established  by 
isogradient  ray  tracing  utilizing  the  five  foot  layer  thicknesses  with  the 
straight  line  segments  as  depicted  (and  exaggerated)  in  Figure  2.  These  rays 
are  used  to  judge  the  rays  formed  by  the  isospeed  method  with  25  foot  layer 
thickness.  The  corresponding  depth-velocity  table  is  formed  by  partitioning 
the  {vi)  into  consecutive  sets  of  size  5  and,  within  each  set,  average  the  five 
values. 

With  the  above  as  background,  the  remainder  of  this  section  deals  with 
numerical  comparisons  treating  three  issues:  the  computational  noise 
generated  by  the  processing  of  a  large  number  of  layers  on  the  computer;  the 
basic  precision  of  the  current  isospeed  ray  tracing;  the  effects  of  extrapolation 
policies  when  the  measured  water  column  does  not  go  sufficiently  deep. 

i.  The  adopted  standard  generally  processes  over  200  water  layers  (20  for 
each  100  feet  separating  source  and  receiver).  The  buildup  of  computational 
noise  can  be  checked  by  using  an  artificial  but  linear  depth  velocity  profile: 
An  exact  ray  can  be  developed  from  the  theory  presented  in  the  previous 
section  and  applied  to  a  single  layer  of  great  thickness.  Then  the  isogradient 
programs  can  attempt  to  match  this  ray  by  tracing  through  the  usual  five  foot 
layers.  This  was  done  for  the  depth  layer  150, 1200  ft.  and  intercepts  4250,  4800, 
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Figure  2.  Schematic  Diagram  Comparing  Isogradient  and  Isospeed 
Representations  of  Depth-Velocity  Information 
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4850  ft/sec  matched  with  slopes  .1,  .05,  .01  ft/secft  respectively.  The  two 
methods  produce  virtually  identical  results.  All  computations  are  in  double 
precision  arithmetic. 

ii.  The  ray  fitting  technique  produces  transit  time,  and  entrance  and  exit 
angles  of  the  ray  connecting  two  points  (source  and  receiver  coordinates)  in 
the  horizontal-depth  plane.  When  isogradient  ray  tracing  is  applied  starting 
at  the  receiver,  using  the  entrance  angle  and  stopping  when  the  transit  term 
is  consumed,  then  the  coordinates  of  the  source  are  reproduced  to  within 
some  small  preassigned  value  e  (we  used  e  =  10    ). 

The  accuracy  of  isospeed  raytracing  using  25  foot  layer  increments  is 
compared  against  isogradient  ray  tracing  using  five  foot-layer  increments.  A 
fixed  set  of  entrance  angle  (60  in  radians)  and  transit  time  (to  in  seconds)  pairs 
have  been  selected.  Generally  they  produce  horizontal  distances  of  1000  to 
6000  feet  and  depths  of  less  than  50  to  over  500  feet.  The  distance,  d,  between 
the  two  versions  of  sound  source  is  compared  for  each  input  pair  (Go,  to)  and 
each  of  the  twelve  DV  profiles.  The  values  are  recorded  in  Table  1  along  with 
the  source  coordinates  (hc,  zc)  for  currently  utilized  isospeed  methodology, 
and  (hs,  zs)  for  the  standard  (isogradient).  An  additional  computation  (hg,  zg) 
was  performed  using  isogradient  ray  tracing  with  25  ft  layer  thicknesses.  (The 
purpose  was  to  provide  an  indication  of  the  relative  importance  of  layer 
thickness  and  the  two  types  of  ray  tracing.)  In  all  cases  the  receiver  depth  is 
1300  feet. 

The  errors  are  computed  using 


=  -\j(hs-hr) 


d  =  \  (hs-hr)2  +  (zs-zz)  (3.2) 
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with  the  subscript  r  replaced  by  c  and  g,  respectively,  in  order  to  identify  the 
current  and  25  ft  isogradient  computations. 

Examination  of  Table  1  shows  that  the  errors,  d,  grow  with  range  and  the 
majority  of  the  total  error  is  in  the  vertical  component.  The  effect  of  varying 
DV  profiles  is  quite  noticeable.  The  larger  errors  can  exceed  one  foot. 
Generally  isogradient  ray  tracing  with  25-foot  layer  thicknesses  is  noticeably 
better  than  the  current  methods. 

iii.  Some  measure  of  the  effect  of  extrapolation  techniques  is  given  in 
Table  2.  The  same  cases  are  developed  as  found  in  Table  1,  but  the  sensor 
depth  has  increased  to  1350  feet.  The  DV  profiles  seldom  goes  below  1300  feet 
so  extrapolation  of  speed  information  is  necessary.  There  are  quite  a  few 
arrays  deeper  than  1300  feet  and  several  (2,  3,  7,  13,  14  see  Table  B-l) 
considerably  so.  Moreover,  the  C-phones  are  even  deeper,  see  Table  B-2. 
(Fourteen  of  the  C-phones  are  deeper  than  1325  ft.)  The  methods  of 
extrapolation  are  explained  in  Appendix  A  along  with  the  visual  effect  of 
their  use.  In  some  instances  the  visual  effect  is  appealing  and  in  others  it  is 
not.  See  the  insets  in  Appendix  A.  Thus  the  values  for  the  standard  (hs,zs)  are 
not  always  well  supported.  Even  so,  the  effect  is  substantial  and  this  is  an 
important  source  of  systematic  error. 

Table  2  is  similar  to  Table  1  in  the  qualitative  sense.  The  effect  of  the  DV 
profile  is  greater  and  at  the  greater  ranges  the  discrepancies  can  be  quite  large. 

4.     ERROR  ASSESSMENT  FOR  THREE  DIMENSIONAL  METHODS 

The  ability  to  locate  a  sound  source  position  from  the  transit  times  (tj,  ..., 
t4)  needs  to  be  assessed  in  three  dimensions  because  of  the  directional 
properties  of  the  array  cubes.   Our  approach  is  to  place  the  acoustic  center  at  a 
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depth  a2,  and  at  the  center  of  a  right  circular  cylinder  of  radius  h.  The  sound 
sources  (k  in  number)  are  equally  spaced  on  a  section  of  the  cylinder  at 
depth  z.  Figure  3  shows  a  plan  view.  Azimuth  is  measured  counter 
clockwise  with  zero  at  "3  o'clock." 

The  ray  fitting  methodology  is  used  to  construct  true  elevation  angles  (81, 
...,  65)  and  true  transit  times  (ti,  ...,  ts)  to  each  of  the  k  sound  sources  from  the 
four  hydrophones  and  the  acoustic  center  (85,  ts).  We  also  need  the  azimuth 
from  the  acoustic  center  (origin  of  circle)  to  the  sources.   These  latter  are 

<J)j  =  27i(j-l)/kforj  =  l/...,k.  (4.1) 

Since  the  ray  fitting  methodology  is  two-dimensional,  we  must 
characterize  the  vertical  planes  connecting  each  sound  source  on  the  cylinder 
to  each  of  the  five  points  in  the  array  cube.  The  technique  for  doing  this  is 
developed  in  Appendix  B.  One  needs  only  the  horizontal  separations  and  the 
vertical  positions  of  source  and  receiver.  Thus  five  rays  are  fit  to  each  sound 
source;  one  to  each  of  the  four  phones  and  one  to  the  acoustic  center.  Also 
five  elevation  angles  are  generated;  but  we  retain  only  the  last,  the  true 
elevation  angle  at  the  acoustic  center. 

During  operations,  the  information  collected  consists  of  the  set  of  transit 
times  (ti,  ...,  t4)  from  sound  source  to  receiver  array,  (see  Figure  1).  For  our 
purposes  these  four  values  are  regarded  as  exact.  Thus  we  can  use  the  values 
produced  in  the  ray  fitting  process.  They  must  be  converted  to  a  ray  tracing 
direction  (azimuth  angle,  (J)  and  elevation  angle,  0)  and  a  single  transit  time 
(tac)  to  stop  the  ray  tracing.  The  currently  employed  procedure  is  described  in 
[5].  It  is  convenient  to  present  certain  aspects  ,  and  this  is  done  in  Appendix  E. 
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Finally  the  values  (0,tac)  together  with  the  array  depth  and  DV  information 
are  fed  into  the  ISOSPEED  ray  tracing  program.  The  output  consists  of  the 
(h,z)  values.  Upon  combining  these  with  the  azimuth  <f>C/  the  estimated 
position  of  the  sound  source  can  be  computed. 

A  number  of  error  calculations  can  be  made.    First  the  error  in  estimating 
transit  time 

timer  =  tac-t5  (4.2) 

affects  the  rule  for  stopping  the  ray  tracing  algorithm.  An  error  of  10"4  seconds 
translates  to  about  half  a  foot  in  horizontal  range,  h,  (i.e.,  u  =  4800  ft. /sec). 
Next,  the  horizontal  error  is  measured  directly 

her  =  hc-hi  (4.3) 

The  values  timer  and  her  should  be  strongly  correlated. 

In  a  similar  fashion,  the  error  in  the  elevation  angle  is  correlated  with  the 
error  in  the  vertical  component 

©er  =  0c  ~  ©o 

(4.4) 
Zer  =  Zc  —  Z\ 

but,  because  of  ray  bending,  the  relationship  is  non-linear.    The  value  zer  is 
measured  directly  in  feet.    The  value  0er  is  more  difficult  to  interpret.   It  is  a 
function  of  the  water  layers  involved.   See  [7]. 
Finally  we  have  the  error  in  azimuth 

tyer  =  <}>c  -  <J>o-  (4.5) 

These  values  are  in  radians  and,  when  multiplied  by  h],  measure  how  far  off 
the  mark  (see  Figure  3)  in  feet,  along  the  cylinder  perimeter. 
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It  is  a  rather  incipient  fact  that  our  errors  are  periodic  functions  of  the  true 
azimuth.  This  is  illustrated  in  Figure  4.  We  suspect  that  this  is  at  the  root 
some  earlier  attempts  to  treat  possible  causes  of  systematic  error,  [6,7]. 


DOWN 
RANGE 


Figure  3.  Cylinder  Cross  Section  Illustrating  Sound  Source 
and  Receiver  Positions 

Many  tabulations  of  these  errors  have  been  made.   Four  sets,  each  with  15 

points  around  the  cylinder  and  four  radii  (hi  =  2000,  3000,  40000,  5000  ft.), 
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Figure  4.  Periodic  Nature  of  Errors  as  a  Function  of  Azimuth 

Case:  DV  of  3/22/89  and  Array  No.  1.  Ref.  Table  3. 
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have  been  selected  for  display  and  they  appear  in  Table  3.  The  basis  of 
selection  was  to  choose  two  of  the  more  interesting  water  columns  (10/27/88 
and  3/22/89,  see  Appendix  A)  each  matched  with  two  of  the  more  severely 
tilted  arrays:  Number  1  and  56,  see  Table  B-l.  The  tilts  for  array  number  1 
have  opposite  signs  while  those  for  array  56  are  of  the  same  sign. 

Let  us  make  a  sample  calculation  in  order  to  aid  in  the  interpretation  of 
these  errors.   The  total  error  offset  is  essentially 

d  =  ^+hl+(hl<t>er)2. 

Thus  for  the  case  of  Array  No.  1  on  10/27/88  at  h}  =  5000  ft.  and  azimuth 
0.418879  radians  (i.e.,  24  deg.  north  of  east)  we  have 

^ (2. 45)2  +  (.  28)2  +  (5  •  2. 3525)2  =12.02  ft.  The  dominant  portion  of  the  error  is 
in  the  azimuth,  which  (measured  along  the  arc)  contributes  2.3525  ft.  of  error 
for  every  1000  ft.  of  horizontal  range.  A  scan  of  Table  3  shows  generally  that 
this  condition  persists,  although  the  periodic  nature  of  the  errors  can  make 
the  effect  small  in  some  directions. 

It  seems  that  the  most  severe  errors  are  associated  with  the  arrays  having 
the  larger  tilts.  We  suspect  that  the  causes  lie  in  the  use  of  the  array  center  as 
the  acoustic  center  and  the  assumption  of  constant  sound  speed  for  the  entire 
array.  The  z-phone  is  about  30  ft  higher  than  the  others.  Also,  the  tilt 
correction  method  is  but  approximate  (see  Appendix  B.) 

We  are  disinclined  to  attach  great  importance  to  the  other  source,  the 
isospeed  ray  tracing.    The  former  was  considered  in  Section  3  and  cannot 
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support  errors  of  this  magnitude.  (The  calculations  in  this  section  apply  the 
polynomial  extrapolation  prior  to  the  development  of  isospeed  layers.  In  this 
way  the  DV  table  extension  is  not  an  issue  in  the  present  comparisons.) 

Accordingly,  the  author  has  proposed  and  tested  a  method  that  treats  the 
two  main  issues.  The  acoustic  center  is  moved  to  the  C-phone.  (Visualize 
Figure  3  with  the  C-phone  on  the  axis  of  the  cylinder.)  In  this  way  the  transit 
time  that  stops  the  ray  tracing  is  t^,  a  directly  measured  quantity.  Also  the 
sound  speed  in  the  layer  containing  the  (X,  Y,  C)  phones  (Z-phone  omitted)  is 
assumed  constant  in  order  to  determine  ((J>p,  0p),  the  azimuth  and  elevation 
angles  of  the  proposed  method.  In  this  way  the  change  in  sound  speed  at  the 
Z-phone  is  taken  out  of  the  computation.   Details  appear  in  Appendix  F. 

Computations  using  this  technique  appear  in  Table  4.  The  improvement 
is  dramatic.  Let  us  repeat  the  earlier  exemplary  computations  10/22/88,  Array 
No.  1,  hT  =  5000  and  <})  =  0.418879.  This  time  the  inputs  come  from  Table  4: 
specifically  V(-22)2  +  (0.5)2  +  (5-0.04306)2  =  0.31  ft.,  a  value  dramatically  smaller 
than  12  ft.  A  scan  of  Table  4  and  comparison  with  Table  3  shows  that  these 
improvements  are  quite  consistent. 

These  computations  also  have  the  advantage  that  the  five-foot  layer 
isogradient  ray  tracing  was  used  together  with  the  exact  tilt  corrections,  eq. 
(B.5.) 

5.     CONCLUSIONS 

Generally  there  are  a  number  of  sources  of  systematic  error.  In  some  the 
effects  are  small,  e.g.,  isogradient  vs.  isospeed  raytracing,  5  foot  versus  25  foot 
layer  thicknesses.  None  the  less  they  are  systematic  and  certainly  no  longer 
necessary.   The  presence  of  systematic  errors  tends  to  build  up  idiosyncracies 
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that  frustrate  the  use  of  standard  statistical  methodology  when 
troubleshooting  other  aspects  of  the  data. 

The  discovery  of  the  periodic  nature  of  the  errors  was  a  surprise.  Its 
presence  adds  another  dimension  to  the  interpretation  of  the  data.  No  doubt 
it  is  the  source  of  some  deception. 

The  major  sources  of  error  related  to  ray  tracing  (see  pages  5,  6)  are 

believed  to  be  associated  with 

(i)     the  conversion  of  transit  times  to  the  direction  of  the  sound 
source 

(ii)    the  constant  speed  extrapolation  of  the  DV  profile  below  1300  ft. 

(iii)  the  use  of  approximations  for  tilt  corrections. 

Of  course  the  effects  of  these  errors  are  directional  because  of  their  periodic 

nature.     We  have  documented  errors  of  12  ft.  or  so  due  to  (i).     We  can 

speculate  10  or  more  feet  because  of  (ii).   Moreover  the  combined  effects  could 

be  additive.  The  effect  of  (iii)  increases  as  the  tilts  increase. 

It  is  further  noted  that  the  above  errors  apply  to  single  arrays.    At  this 

point  we  have  no  comment  about  how   they  may  combine  to  produce 

mismatches  in  the  array  overlap  regions.    It  would  be  more  comfortable  to 

treat  this  issue  by  introducing  the  changes  and  then  collecting  more  data. 


6.     ADDENDUM 

It  came  to  the  author's  attention  during  the  final  editing  phase  of  this 
report,  that  the  approximate  tilt  corrections,  eq.  (B.4),  are  not  the  ones 
currently  employed  by  the  software.   Rather,  the  system  is  using 


(6.1) 


X(l)' 

c2      0      _s2 

'Xo(i)' 

X(2) 

>  =  < 

0     q     -s-[  ' 

'  X0(2) 

X(3). 

S2      Sj      CyC2 

*o(3). 
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Accordingly    the    author    made    some    specialized    computations    for 

purposes  of  indicating  the  effect.   Referring  to  Tables  3  and  4,  we  have  chosen 

to  treat  the  exemplary  case:   water  column  of  10/27/88,  array  number  1,  \\\  = 

5000  ft,  zi  =  250,  and  $  =  0.418879  radians.   The  table  below  contains  the  four 

errors  and  total  offset  for  each  of  four  algorithms: 
(i)     the  isospeed  method  (from  Table  3) 
(ii)    the  isospeed  method  using  (6.1)  vice  (B.4) 
(iii)  the  isospeed  method  using  (B.5)  vice  (B.4) 
(iv)  the  proposed  method  (from  Table  4) 

TABLE  5.  EFFECT  OF  TILT  CORRECTION  METHOD 

10/27/88,  Array  #1,  hj  =  5000  ft,  Zi  =  250  ft,  (J)  =  0.418879  rad.,  a2  =  1308.76  ft. 


<j)er  Qer  ner 


-er 


(i) 

.0023525 

.0003869 

.28 

-2.45 

12.02 

(ii) 

.0021174 

-.0001567 

1.47 

2.95 

11.09 

(iii) 

-.0000161 

-.0000077 

1.26 

2.19 

2.53 

(iv) 

.0000166 

.0000430 

-.05 

-.22 

0.31 

The  results  of  Table  5  suggest  the  following:  there  is  little  distinction 
between  the  use  of  (B.4)  and  (6.1)  for  the  tilt  correcting  (compare  (i)  and  (ii)). 
However  case  (iii)  suggests  that  there  is  much  to  be  gained  by  use  of  the  exact 
tilt  correct  (B.5)  even  if  isospeed  ray  tracing  is  used  with  25  ft  layer  thickness. 
Finally  case  (iv)  suggests  that  considerable  gains  are  available  if,  in  addition, 
we  use  the  proposed  methodology. 

All  of  these  systematic  errors  are  mathematical,  i.e.,  due  to  choice  of 
algorithms.   There  is  no  longer  any  reason  not  to  use  the  best. 
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APPENDIX  A 

The  twelve  depth  velocity  profiles  used  in  this  study  are  recorded  here  in 
graphical  form.  They  include  the  profiles  used  in  [7].  Generally  the  measured 
values  stop  at  a  depth  of  about  1300  feet  and  it  is  necessary  to  extrapolate,  in 
several  instances  to  depths  greater  than  1350  feet  (see  Table  B-2). 

The  insets  of  the  graphs  illustrate  two  different  extrapolation  schemes. 
The  constant  value  extrapolation  is  the  one  currently  in  use,  [5].  But  there  is  a 
slight  deception.  Current  methodology  utilizes  isospeed  profiles,  see  Figure  2, 
and  25  foot  layer  thicknesses.  The  constant  value  extrapolation  is  the  value  of 
the  deepest  25  foot  layer  appearing  in  an  isospeed  profile. 

The  other  extrapolation  is  based  upon  fitting  a  second  order  polynomial 
to  the  deepest  hundred  feet  of  the  original  profile  (five  foot  layer  thicknesses). 
There  are  two  steps  in  this  process.  First  is  fitting  the  curve  by  least  squares. 
Because  of  the  equally  spaced  depth  increments,  the  fitting  takes  an  especially 
simple  form.    Using  the  equation, 

v  =  a  +  b(u-u)  +  c(u-u)  (A.l) 

where  u  is  velocity,  u  is  layer  depth,  and  u   -  average  depth,  the  normal 
equations  take  the  form 

Zv  =  Za -  Zb(u  -  u) -  cZ(u  -  u) 
Zv(u-u)  =  Za(u-u)- Zb{u-u)   -cZ(u-Tiy 
Zv{u  - u)2  =  Za(u  - u)2  -  Zb{u  -  u )3  - cZ(u  -  u)4 


37 


Sv 

=  na  +  cSu2 

=  bSu2 

s 

=  aSu2  +  cSt 

Using  the  notation  Sv  =  £v,  Suu  =  £v  (u-  u  ),  Suuu  =  Lv  (u-  u  )2.  SU2  =  Z  (u- 
u  )2,  SU4  =  Z(u-  u  )4and  recognizing  Z(u-u)=  I(u-u)    =0  because  of  the 
uniform  spacing,  the  above  equations  assume  the  reduced  form 

(A.2) 

and  we  solve  for  the  coefficient  of  the  quadratic  term 

c  =  inSvuu  -  SvSu2)/  (nSuA  ~  ^2)  (A3) 

(This  done,  values  for  a  and  b  come  easily  from  the  first  two  equations.) 

We  note  in  passing  that  if  n  =  2K+1,  an  odd  number,  and  A  is  the  spacing 
between  consecutive  values  of  {u},  then  u  is  an  integer  and 

A2 

Sul  =  \[m^)f,  (A.4) 

nSw4  "  Su2  =  ^4K2(K  +  1)2{2K  + 1)(7  -  4K)  /  18 

Such  formulas  expedite  the  computation  of  c  in  (A. 3). 

The  second  step  deals  with  the  issue  of  how  to  use  this  quadratic  for 
extrapolation  purposes.  Generally,  the  direct  use  of  the  equation  (A.l)  would 
lead  to  a  visual  discontinuity  between  the  last  measured  value  and  the  first 
extrapolated  one.  We  have  chosen  not  to  do  this.  Instead  we  recognize  that 
Ac  represents  the  first  difference  of  the  gradient  sequence  {t>i(i)},  see  Figure  2. 
So  to  preserve  continuity,  the  extrapolation  is  enabled  by  using  successive 
updates 
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U](z)=  n[(i-l)  +  A  -c 
v(i  +  l)-  v(i)  +  Av\(i) 
vQ(i)  =  [u(i  +  l)v(i)-u(i)v(i+l)]  I A 

where  \)(i)  is  the  estimated  sound  speed  at  u(i);  \>o(i)  and  \)i(i)  are  the  intercept 
and  slope  values  of  the  straight  line  fit  for  the  ith  layer.  It  is  this  continuity 
preserving  step  that  explains  the  occasional  appearance  of  "crooked" 
extrapolations  in  the  insets. 

The  use  of  isogradient  ray  tracing  with  25  foot  layer  thicknesses  was  used 
in  Section  3.  If  extrapolation  of  the  water  column  was  required,  then  a 
different  second  order  polynomial  method  was  used,  and  after  the  conversion 
to  25  foot  layers.  Basically  the  quantity  Ac  was  estimated  by  averaging  the 
difference  of  the  last  five  values  of  x>\  (last  100  feet),  and  then  proceeding  in 
the  same  way  as  stated  above.  No  graphs  showing  the  effect  of  this  have  been 
prepared.  But  the  choice  does  have  an  effect  upon  the  dg  values  computed  for 
Table  2. 

We  take  this  opportunity  to  note  that  the  visual  appeal  of  the  quadratic 
extrapolation  is  rather  good  in  instances  1,  2,  3,  4,  6  and  11.  The  others  are 
easily  challenged.  This  information  may  influence  the  reader's  interpretation 
of  some  values  appearing  in  Table  2. 
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APPENDIX  B 


COORDINATE  SYSTEMS 

The  developments  that  follow  deal  with  several  coordinate  systems  and 
we  need  efficient  ways  to  distinguish  among  them.  First  of  all  is  the  right 
handed  system  defined  by  the  X,  Y,  Z  and  C  hydrophones  of  the  array.  All 
incoming  transit  times  must  be  interpreted  in  this  system  and  we  call  it  cs(a), 
or  the  coordinate  system  of  the  array.  The  origin  is  at  the  c  hydrophone. 

Since  this  system  is  generally  tilted  with  respect  to  a  "flat  earth  surface" 
system  there  is  need  to  rotate  cs(a)  into  alignment  with  a  common  coordinate 
system  for  all  arrays,  or  a  range  coordinate  which  has  horizontal  directions 
consistent  with  the  earth's  surface.  Such  a  resultant  system  will  be  called  cs(b) 
and,  for  convenience  of  terminology,  the  X  arm  of  cs(a)  is  rotated  to  a  position 
called  east;  the  Y  arm  to  a  position  called  north,  and  the  Z  arm  to  a  position 
called  vertical  or  zenith.   The  origin  is  still  at  the  c  hydrophone. 

The  ray  tracing  methodology  of  Ref  [5]  attempts  to  locate  the  sound  source 
in  relation  to  a  specific  point  termed  the  acoustic  center.  This  center  is  the 
geometric  center  of  the  array  cube.  The  resulting  coordinate  system  is  a 
translation  of  cs(a),  and  will  be  called  cs(ac).  In  this  system  the  four 
hydrophone  positions  are  specified  by  D/2  times  the  vectors 
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respectively,  whereas  in  cs(a)  these  four  vectors  would  be  the  three  unit  basis 
vectors  and  the  origin.  Conversion  from  cs(ac)  to  cs(a)  is  a  direct  translation. 

The  conversion  of  a  point  located  in  cs(a)  to  the  cs(b)  system  requires  a 
three-dimensional  rotation  based  upon  the  tilt  angles  and  the  "ZROT" 
horizontal  direction  correction.  It  is  convenient  to  describe  this  in  terms  of 
the  three  Euler  angles  <{>i,  (j>2,  <J>3,  or  roll,  pitch  and  yaw.   Letting 

sj  =  sine  (<j>j)     and     Cj  =  cosine  (<)>i)       i  =  1,  2, 3  (B.2) 

we  define  three  successive  rotations  which,  when  applied  sequentially  to  a 
point  in  cs(a),  will  in  the  end  describe  it  in  cs(b).  First  hold  the  X  arm  fixed 
and  rotate  the  Y-Z  plane  through  an  angle  <j>i;  the  matrix  of  this  transform  is 

1     0      0 

0    ci    -Si 

L  0    si      ci    J 


Pi 


Next  hold  the  (current  position  of)  Y-arm  fixed  and  rotate  the  X-Z  plane 
through  an  angle  of  $2;  this  transformation  has  matrix 

C2    0    -S2 
0     1      0 

_   S2      0       C2     _ 


P2  = 


Finally  hold  the  (current  position  of)  Z-arm  fixed  and  rotate  the  X-Y  plane 
through  an  angle  of  <j>3;  the  transform  is 

C3      -S3      0 

p3=     s3      c3     0 
L   0       0      1  J 

The   successive   applications   of   these   three   rotations   is   a   unitary 
transformation,  (i.e.,  its  inverse  is  equal  to  its  transpose), 
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B  =  p3  P2P1  (B.3) 

and  if  a  is  a  three  vector  in  cs(a),  then  b  =  Ba  is  the  same  vector  referenced  in 
cs(b). 

The  determination  of  the  Euler  angles  is  accomplished  as  follows.  The 
submerged  arrays  have  tilt  indicators  on  the  X  and  Y  arms  which, 
individually,  measure  the  angles  that  these  arms  make  with  the  horizontal. 
An  accounting  for  these  tilts  must  be  made  when  the  ray  trace  azimuth  and 
elevation  angles  are  converted  to  a  horizontal  based  coordinate  system.  The 
technique  currently  in  use  takes  the  apparent  position,  Xo,  and  applies  the 
transformation 


(B.4) 


so  that  the  new  apparent  position,  X,  is  referenced  in  a  plane  level  with  the 
earth.  This  transformation  is  an  approximation  which  simply  rotates  the  two 
arms  to  the  horizontal  as  if  they  were  separate  unconnected  arms  and  the 
rotation  of  one  does  not  affect  the  rotation  of  the  other.  That  is,  the  first  two 
rows  of  the  coefficient  matrix  are  not  orthogonal.  The  result  is  an 
approximation  whose  success  depends  upon  the  smallness  of  the  tilt  angles. 

The  exact  way  to  accomplish  this  goal  involves  the  direct  replacement  of 
the  coefficient  matrix  with  the  product  rotation  p2pi: 
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Upon  comparing  these  two  coefficient  matrices,  one  sees  there  is  choice  in 
identifying  the  two  tilt  angles  with  these  two  Euler  angles.  We  have  chosen 
to  match  the  first  two  elements  of  the  third  row: 

s2  =  sin(XTILT)        sA  =  sin(YTILT)  /  c2.  (B.6) 

The  geometric  interpretation  is  as  follows.  First  hold  the  X  arm  fixed  and 
rotate  the  plane  of  the  Y-Z  arms  so  that  the  Y  arm  is  horizontal.  This  is  not  a 
vertical  projection.  The  division  by  C2  shows  that  one  must  rotate  through  an 
angle  greater  than  YTILT  in  order  to  maintain  orthogonality  of  the  coordinate 
system  when  making  the  Y  arm  parallel  to  the  earth's  surface.  This  done,  we 
next  hold  the  Y  arm  fixed  in  its  new  position  and  rotate  the  plane  of  the  X-Z 
arms  so  that  the  X  arm  is  horizontal.  Since  the  new  Y  arm  is  already 
horizontal  this  second  rotation  is  a  vertical  projection  through  an  angle  of 
XTILT. 

This  latter  method  is  exact.  The  nature  of  the  original  approximation  can 
be  assessed  by  comparing  the  two  coefficient  matrices,  using  numerical 
inputs.  The  effect  is  not  great  for  most  of  the  tilts  present  at  Nanoose. 

To  complete  the  conversion  of  cs(a)  to  cs(b)  we  must  find  the  third  Euler 
angle  in  terms  of  ZROT,  the  rotation  of  the  horizontal  plane  determined  by 
array  survey.  In  [5],  ZROT  is  defined  as  the  angle  from  the  horizontally 
projected  X  arm  to  the  range  center  line  (east).  The  comparison  of  this 
definition  with  that  of  03  (see  p3>  leads  to 

s3  =  -sin(ZROT)  (B.7) 

The  subsequent  application  of  p3  to  P2P1  will  bring  the  point  a  into  east-north 
orientation. 
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Finally  the  position  location  system  prefers  to  specify  an  object's  position 
in  terms  of  east,  north,  and  depth  (positive)  below  the  sea  surface.  This 
system  is  a  left-handed  one  with  origin  on  the  sea  surface  directly  above  the 
acoustic  center  of  the  array. 

Table  B-l  contains  the  positions  of  the  Nanoose  arrays  at  the  time  of  their 
most  recent  survey  dates.  These  are  the  positions  of  the  acoustic  centers  in 
the  range  coordinate  system.  It  is  both  useful  and  instructive  to  use  our 
coordinate  system  superstructure  and  locate  the  hydrophones  of  each  array  in 
the  range  coordinate  system. 

Let  a  be  the  3-vector  locating  the  acoustic  center  of  an  array  in  the  range 
coordinate  system,  see  Table  B-l.  Let  a  be  the  location  of  one  of  the  phones  in 
cs(a)  and  let  f, 


f=(D/2) 


1 


(B.8) 


be  the  position  of  the  acoustic  center  in  cs(a).  In  the  system  cs(b)  these  vectors 
are  Ba  and  Bf,  and  the  location  of  the  phone  relative  to  the  acoustic  center  is 
B(a-f).   The  position  p  of  the  phone  in  the  range  coordinate  system  is 

p  =  oc+B(a-f)  (B.9) 

Since  the  C-phone  is  the  origin  in  cs(a),  its  range  coordinate  position  is  a  -  Bf. 
the  result  of  this  computation  is  in  Table  B-2. 
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TABLE  B-l.  NANOOSE  ARRAY  COORDINATES  AND  ORIENTATION 

ANGLES 

(Distances  in  feet;  angles  in  radians) 

X 

12188.01 
19463.16 
26991.39 


Survey  Date 

AR  0  6/20/85 
AR 1  6/20/85 
AR  2         7/12/85 


AR3  1/7/88  34505.10 

AR4  10/24/88  42005.19 

AR  5  6/20/85  49497.00 

AR6  6/20/85  56972.28 

AR  7  7/30/85  64680.66 

AR8  11/16/88  71969.73 

AR  9  5/7/84  3.00 

AR  10  3/12/84  47100.00 

AR11  7/18/85  23173.89 

AR12  6/20/85  30731.25 

AR  13  6/20/85  38213.61 

AR  14  6/20/85  45647.07 

AR  15  6/19/85  53249.43 

AR  16  9/13/85  60859.74 

AR  17  6/16/87  68217.93 

AR  54  2/2/88  38029.95 

AR  55  6/20/85  45645.75 

AR56  7/30/85  53180.13 

AR57  7/30/85  60745.71 

AR23  6/20/85  41605.14 

AR  24  4/17/89  49572.00 

AR  25  10/24/88  56993.79 

AR  26  8/8/88  64442.94 

AR27  7/15/80  22119.60 

AR  28  5/4/83  45000.00 

AR  29  2/2/79  0.00 


Y 

Z 

XTILT 

YTILT 

ZROT 

-131.52 

-1295.33 

0.002909 

0.014835 

-0.208183 

-174.99 

-1308.76 

0.061523 

-0.036070 

1.362579 

-109.83 

-1323.25 

0.000145 

0.005236 

2.670336 

-80.76 

-1323.32 

0.027925 

-0.011345 

2.928139 

-55.17 

-1318.28 

0.001164 

-0.040288 

-2.315877 

-25.23 

-1315.58 

-0.000291 

-0.004072 

1.668535 

-21.21 

-1308.50 

0.013817 

0.041161 

-0.703420 

15.33 

-1353.39 

0.034907 

0.022835 

-0.574144 

-29.28 

-1300.89 

-0.005963 

-0.012217 

-1.577341 

3.00 

1.00 

0.000000 

0.000000 

0.000000 

-3600.00 

-1300.00 

0.000000 

0.000000 

0.000000 

-6488.40 

-1312.09 

-€.004654 

0.000436 

2.784376 

-6553.05 

-1312.90 

0.002036 

0.001745 

-3.042179 

-6640.77 

-1323.05 

0.000291 

0.006254 

1.373522 

-6513.18 

-1324.78 

0.001309 

0.002327 

-2.348044 

-6354.60 

-1316.66 

0.003345 

0.004509 

0.581544 

-6356.07 

-1313.42 

0.014835 

0.036943 

2.303276 

-6524.10 

-1313.43 

0.008290 

0.034761 

2.158449 

5401.98 

-1212.69 

0.007709 

-0.003782 

-1.056919 

6369.66 

-1188.12 

0.027634 

0.039415 

-0.728553 

6417.96 

-1218.84 

0.037525 

0.048142 

-1.392651 

6419.40 

-1088.24 

0.006981 

0.001891 

-3.108606 

12150.18 

-1268.23 

-.002182 

0.003200 

-1.845214 

-12966.00 

-1300.00 

-0.007272 

0.055269 

-1.343904 

-12999.33 

-1205.48 

0.000291 

-0.002182 

-0.593726 

-12971.04 

-1255.35 

-0.014835 

-0.012654 

3.134192 

-15908.70 

83.00 

0.000000 

0.000000 

0.000000 

1500.00 

-1350.00 

0.000000 

0.000000 

0.000000 

0.00 

0.00 

0.000000 

0.000000 

0.000000 
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TABLE  B-2.  LOCATIONS  OF  THE  C-HYDROPHONES 


AR 

Date 

XC 

YC 

zc 

0 

6/20/85 

12170.190 

-143.316 

-1310.105 

1 

6/20/85 

19473.791 

-192.189 

-1325.075 

2 

7/12/85 

27011.563 

-103.354 

-1338.287 

3 

1/7/88 

34522.511 

-69.103 

-1338.681 

4 

10/24/88 

42004.318 

-33.398 

-1332.413 

5 

6/20/85 

49513.397 

-38.633 

-1330.630 

6 

6/20/85 

56950.933 

-23.551 

-1323.123 

7 

7/30/85 

64659.408 

10.557 

-1367.552 

8 

11/16/88 

71954.918 

-13.999 

-1315.793 

9 

5/7/84 

-12.000 

-12.000 

-14.000 

10 

3/12/84 

47085.000 

-3615.000 

-1315.000 

11 

7/18/85 

23193.258 

-6479.598 

-1327.004 

12 

6/20/85 

30744.657 

^6536.662 

-1327.956 

13 

6/20/85 

38225.375 

-6658.513 

-1337.943 

14 

6/20/85 

45646.877 

-6492.002 

-1339.829 

15 

6/19/85 

53245.085 

-6375.441 

-1331.552 

16 

9/13/85 

60880.699 

-6357.757 

-1328.681 

17 

6/16/87 

68238.605 

^6528.792 

-1328.448 

54 

2/2/88 

38009.399 

5407.725 

-1227.510 

55 

6/20/85 

45624.165 

6367.887 

-1202.470 

56 

7/30/85 

53162.159 

6429.360 

-1233.742 

57 

7/30/85 

60760.102 

6434.858 

-1103.370 

23 

3/20/85 

41594.799 

-12131.725 

-1283.312 

24 

4/17/89 

49554.120 

-12955.612 

-1315.728 

25 

10/24/88 

56972.961 

-13003.338 

-1220.483 

26 

8/8/88 

64458.271 

-12955.966 

-1269.935 

27 

7/15/80 

22108.374 

-15890.700 

68.000 

28 

5/4/83 

44985.000 

1485.000 

-1365.000 

29 

2/2/79 

-15.000 

-15.000 

-15.000 
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APPENDIX  C 

2.     RAY  FITTING  AND  RAY  TRACING 

A  "ping"  sound  source  produces  a  wave  front  that  travels  through  the 
water  and  is  detectable  by  the  receiving  transducers.  A  ray  is  the  path 
generated  by  the  normal  to  the  wave  front.  Ray  tracing  is  the  activity  of 
following  a  ray  from  a  receiver  for  a  fixed  amount  of  time  in  order  to  locate 
the  sound  source.  Ray  fitting  is  the  activity  of  recreating  the  ray  path  from  the 
positions  of  the  source  and  the  receiving  sensor.  The  latter  is  needed  to  study 
the  error  performance  of  the  former. 

The  speed  of  sound  in  water  is  assumed  to  vary  with  depth,  but  remain 
homogeneous  in  the  horizontal.  Thus  a  single  depth-velocity  profile  is  valid 
throughout  the  field. 

We  begin  with  some  preliminaries  followed  by  the  development  of  a  ray 
fitting  algorithm,  which  may  be  viewed  as  an  inverse  to  ray  tracing.  It  is 
certainly  more  difficult.  All  paths  are  direct  paths.  That  is,  no  provision  is 
made  for  reflections  or  refractions  that  produce  non  monotone  rays 
connecting  source  and  receiver.  Also  our  interest  lies  in  the  greater  ranges 
and  no  adjustments  are  presented  for  sources  that  are  directly  above  the 
receiver. 

The  notation  is  chosen  to  be  consistent  with  that  used  in  the  Fortran 

source  codes.   We  are  given  a  speed-depth  profile  in  the  form  of  pairs  lm,  and 

velj. 

lmj  =  depth  of  ith  water  layer,  positive  down, 
velj  =  speed  of  sound  at  the  depth  lm,. 
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Digression.  The  water  speed  processing  system  at  NUWES  produces 
average  velocity  values  for  a  large  number  of  equally  spaced  points.  These 
averages  are  intended  to  serve  as  constant  values  for  the  entire  layer,  eq.  (3.1). 
Thus  the  information  is  provided  in  the  form,  for  i  =  1,  ...,  m 

li  =  depth  of  water  layer  boundary 
velj  =  velocity  constant  for  layer  (lj,  lj+i) 

The  values  used  for  lmj  above  are  the  layer  midpoints: 
lmj  =  .5  x  (lj  +  lj+i)  for  i=l, ...,  m-1 

The  algorithms  developed  below  do  not  require  layers  of  equal  thickness. 
Thus  they  can  accommodate  the  user  who  wants  to  use  thin  layers  at  depths 
of  rapid  change  and  thick  layers  at  depths  of  slow  changes.  All  computations 
should  be  in  double  precision  arithmetic. 

Occasionally  as  application  leads  to  a  receiver  whose  depth  is  greater  than 
the  largest  {lmj}.  For  these  cases  we  have  been  using  an  extrapolation  of  the 
sound  speed  profile  that  adjoins  a  sufficient  number  of  layers,  all  of  the  same 
thickness  as  the  deepest  of  the  original  layers.  The  corresponding  velocities 
are  extrapolated  using  a  second  order  polynomial  calculated  from  a  fixed 
second  difference  whose  value  is  the  average  of  the  four  deepest  second 
differences.  (Use  of  the  coefficient  of  the  quadratic  term  in  a  least  square  fit 
has  also  been  employed  as  an  option.) 

Since  ray  fitting  is  a  rather  delicate  operation  we  use  the  isogradient 
technique.  I.e.,  straight  lines  are  fitted  to  the  speed  for  each  layer;  the  constant 
gradient  of  slope  is  computed;  the  profile  itself  is  then  a  continuous  function 
of  connected  straight  line  segments.   See  Figure  2.   I.e.,  if 

lmj  <  z  <  lirti+i 
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then 

vel(z)  =  vo(i)  +  vi(i)z 

where 

vo(i)  =  lmi+i  •  veli-lmj  •  veli+i 

vi(i)  =  (veli+i-veli)/dzi 

dzj  =  lm  i+i-lmj 

Typically  the  values  of  (vo(i)}  are  positive  and  large  compared  to  the 
(vi(i)}  which  are  small  and  of  either  sign  or  zero.  Snell's  law  comes  into  play 
and  the  ray  invariant  will  be  denoted  rv.  The  ray  path  in  each  layer  will  be  a 
circle  arc  (vi(i)^O)  or  a  straight  line  (vi(i)=0). 

Now  we  are  positioned  to  describe  our  ray  fitting  algorithm.  It  is  two 
dimensional  (horizontal,  vertical)  and  given  the  endpoints 

(ai,a2>    receiver 
(pi/P2)    sound  source 

the  goal  is  to  compute 

00  entrance  angle  at  the  receiver 

01  exit  angle  at  the  sound  source 

t        transit  time  of  sound  from  source  to  receiver 

There  is  no  loss  in  placing  the  origin  on  the  water  surface  directly  over  the 
receiver.  Thus,  ai  =  0  and  p2  <  &2  (depth  of  receiver). 

The  algorithm  is  an  iterative  one  that  operates  as  follows.  Initialize  the 
entrance  angle  0o-  Use  this  angle  to  ray  trace  upward  through  the  layers  until 
the  vertical  value  of  p2  is  achieved.    Compute  the  current  horizontal  value  h, 
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and  compare  it  with  pi.    If  this  is  within  a  preassigned  small  number,  e,  stop 
Figure  C-l.  Details  of  single  layer  processing 

v(z)  =  vo  +  viz     (linear  profile)     Zi<z<zo 

Given:  60,  zo,  zT/  v0/  Vi,  and 

•  cosQo) 

the  ray  invariant  rv  =  — ; — r~ 
3  v(z0) 

Case  vi  =  0 

dz  =  z\  -  zq 

dw  =  dz/sin(Oo) 

hi  =  ho  +  dw  •  cos(0o) 

dt  =  dw/vo 

cos(6i)  =  rv  •  v(zi) 

Case  vi  *  0 
q2  =  -vo/vi 
s  =  sin(6o) 
c  =  cos(0o) 
qi  =  ho  +  (q2-z0)s/c 
r  =  signum(q2)(q2  -  z0)/c 
h  =  q!  -  signum(q2)  r  •  sin(9) 

1  ffli    dd 


dt  =  ±\ 

7»      J( 


v-i  J°o  cos(0) 


=  lln 


^Ucin/fl\^ 


*>1 


v 


l  +  sin(fl) 
cos(0) 
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CO;>(8i)  =  rv  •  v(Z]) 

and  compute  the  transit  time.     Otherwise,  adjust  the  entrance  angle  60 

according  to  the  ratio  of  the  current  rise  over  run, ,  to  the  desired  one, 

a2-p2 

Pi 

,  and  repeat  the  algorithm  using  the  new  initialization. 

ALGORITHM  RAYFIT 

Initialize 

(i)  Determine  the  layers  that  contain  the  source  and  receiver. 
Choose  j,  n  so  that 

lmj  <  p2  <  lmj+i 

lmn  <  a2  <  lmn+1 
(ii)    Make  thickness  corrections  for  the  extreme  layers 

dzj  =  lmj+!-p2 

dzn  =  a2~lmn 
(iii)  Compute  the  sound  speed  at  depths  a^  p2 

va2  =  vo(n)  +  vi(n)  •  a2 

vp2  =  vo(j)  +  vt(j)  •  p2 

(iv)  Unless  a  "previous"  value  for  60  is  available,  fit  a  straight  line 
through  the  depth-velocity  profile  in  the  range  (p2,  a2)  and  use 
Go  corresponding  to  the  circle  arc  (or  line)  of  that  approximate 
profile.  I.e., 

If  va2  =  vp2  then  60  =  tan_1((a2  =  P2VP1),  else 

va2  *  P2-VP2  '  a2 

q2=_      (a2-P2) 
qi  =  .5  •  pi  +  .5  ■  (p2-a2)(p2+a2-2q2)/pi 
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6o  =  tan-1  {qi/(q2-a2)l 
endif 
Set  initial  values  for  iteration 

A.  i  =  n,  s  =  sin(0o),  c  =  Vl-s 
rv  =  c/va^                   h  =  0,                           z  =  a2 

Main  raytracing  code 

B.  If  vi(i)  =  0,  then 
dw  =  dzj/s 

h  =  h  +  c  • dw 
else 

q2  =  -vo(i)/vi(i) 

qi  =h  +  (q2-z)-  s/c 

r  =  V(h-qi)2  +  (z-q2)2 
c  =  rv  •  vel(i) 

s  =  VT:c2 

h  =  qi  -  signum(q2)  •  r  •  s 
endif 

If  i  =  j,  goto  TEST 

z  =  lmj 

i  =i-l 

goto  B 
TEST:    0i  =  cos_1(rv  •  vp2) 

If  V](j)  *  0         h=qi-signum(q2)  •  r  •  sin(Oi) 

If  I  h-p2 1  <  e ,  goto  FINI 

Re-estimate  9o 

90  =  tan-MtanOo)  •  h/pi) 

goto  A 
FINI: 

ang(j)  =  Gi 

angn+1=80 

angi  =  cos_1(rv  •  vel(i))  for  i  =  j+1,  ...,  n 
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Compute  transit  time 

t  =  0 

Do  TC    i  =  j,  n 

If  vi(i)  =  0 

t  =  t  +  dzj/(v0(i)  •  sin(angi)) 

else 

1        Jcos(angi+i)(l-fsin(angi)) 
* "  vi(i)   n[(l+sin(angi+1))cos(angi) 

endif 

TC    continue 

Remove  extreme  layer  thickness  corrections 
dzi  =  lnrtj+i  -  lmj 
dzn  =  lmn+i  -  lmn 
end 

This  algorithm  has  been  quite  useful  to  us.  Using  e  =  10-6  we  typically 
have  8  to  10  iterations  through  A. 

Raytracing  algorithms  are  less  sensitive,  but  since  a  good  one  is  readily 
available  by  merely  modifying  the  above,  let  us  do  that.  The  process  is  an 
inverse  one  in  that  the  goal  is  to  compute  pi,  p2  given  6o  and  to,  where  to  is 
the  transit  time. 

ISOGRAD 

Initialize  by  locating  the  layer  containing  the  receiver,  establishing  the  ray 
invariant,  etc. 

Choose  n  so  that  lmn  <  a2  <lmn+i 

i  =  n,  h  =  0,  z  =  a2,  t  =  0 

s  =  sin(Go),  c  =  cosOo) 

va2  =  v0(i)  +  v^i)  •  a2 
rv  =  c/va2 
dzn  =  a2  -  lmn 
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AA:    If  vi(i)  =  0,then 

dw  =  dzj/s 
dt  =  dw/vo(i) 
h  =  h  +  c • dw 

else 

q2  =  -vo(0/v1(i) 
qi  =  h  +  (q2-z)s/c 

r  =  V(h-qi)2  +  (z-q2)2 
cp  =  rv/vi(i) 

sp  =  Vl-cp2 

l       r  c  l+sp 
dt  =  ^i)lniT^  cp 
h  =  qi-signum(q2)  sp 

endif 

t  =  t  +  dt 

If  t  >  t0  goto  FINAL 
z  =  lrrtj 

s  =  sp  c=cp  i  =  i-1  goto  AA 

FINAL 

dt  =  to  +  dt  - 1 
lfvi(i)  =  0 
dzj=  vo(i)  •  dt 
dw  =  dzj/sp 
pi  =  h  +  dw  •  cs 
P2  =  z-  dZj 
else 

x  =  exp{dt-  vi(i))(l+s)/c 
cp  =  2x/(l+x2) 

sp  =  Vl-cp2 
pi  =  qi  +  r  •  sp 
p2  =  q2  +  r  •  cp 
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endif 

Restore  layer  integrity 

dzj  =  lrrtj+i  -  lrrii 
dzn  =  lmn+i  -  lmn 

Compute  exit  angle 
01  =  cos_1(cp) 
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APPENDIX  D.  EFFECT  OF  SOUND  SPEED  OSCILLATION  WITHIN  A 

LAYER  OF  WATER 

Suppose  a  water  layer  of  thickness  Az  has  an  average  sound  speed  of  \)  feet 
per  second.  Suppose  further  that  the  speed  profile  within  the  layer  is  an 
oscillation  of  frequency  K  and  amplitude  5.  We  address  the  question  of  how 
this  profile  can  affect  the  nominal  calculation  of  the  ray's  horizontal  distance 
and  transit  time  through  the  layer. 

The  question  is  most  easily  treated  using  isospeed  ray  tracing  and 
modeling  the  oscillations  as  follows:  Partition  the  layer  into  2K  equithick 
sublayers  and  assume  the  sound  speed  alternates  between  the  values  u+8  and 
■0-8  as  we  move  through  these  layers.  Let  0i  be  the  elevation  angle  for  the 
"0+5  speed  layers  and  02  for  the  \)-8  layers. 

According  to  isospeed  ray  tracing  formulation,  the  horizontal  distance 
advanced  in  the  layer  is 

2K 


and  the  transit  time 

1 1/  *—t 


A>«        1 


ij;-sin(0;-) 


Az 


I 


(v  +  8)sm(6^)     {v-8)sm{62) 


The  nominal  values  are  found  by  using  \)  as  the  speed  and  0i  as  the  angle 
throughout  the  layer.    Thus  the  error  in  these  two  values  is 
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ah  =  ^[cot(e2)-cot(e,)]f 


AT  =  ^ 


2  [(v+S)sm(0i)     (v-8)sin(d2)     usin^)] 

and  9i  and  02  are  related  by  the  ray  invariant  equation 

cos(0j)  _  cos(#2) 
v  +  8         v-S 


The  error  AH  is  affected  by  velocity  only  through  this  equation.  Notice  that 
the  errors  do  not  depend  upon  the  frequency  K.  Using  8/u  as  the  proportion 
of  the  speed  appearing  in  the  amplitude,  we  can  rewrite 


4T  =  ^ 
2v 


+  *z 


+ 


2v 


2v 


(l  +  8/v)sin(9^)     (l-<5/u)sin(02)     sin(^) 

1-8/v     1+8/v         2 

: ^ : 

sinffy)       sin(02)      sinffy) 

f  \ 


1 


1 


sin(#2)     sinffy) 


(1  +  8/  v) 


since  the  proportion  8/u  is  believed  small. 

Some  speculative  calculations  appear  in  Table  (D-l)  for  Az  =  5  feet  and  u  = 
4800  feet/ second. 

The  calculation  is  not  very  sensitive  to  values  of  \),  but  quite  responsive 
to  the  elevation  angle.  It  should  be  noted  that  the  signs  of  AH,  AT  change  if 
the  modeled  oscillations  are  in  reverse  order.  This  is  equivalent  to  replacing 
with  -5. 

This  author  is  not  qualified  to  judge  the  reality  of  the  suggested  values  of 
5/u   Certainly  the  question  deserves  more  attention.   It  is  common  to  process 
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some  200  of  these  5  foot  layers  in  a  ray  computation.   The  error  buildup,  even 
if  the  signs  change  randomly,  could  be  significant,  perhaps  15  times  AH. 

TABLE  D-l.  EFFECT  OF  OSCILLATION 

(v  =  4800  feet/second;  Az  =  5  feet) 


5A) 

6i 

AH 

AT 

.0001 

0.1 

0.48 

-.0001006 

0.15 

0.15 

-.0000301 

0.2 

0.06 

-.0000127 

0.3 

0.02 

-.0000037 

0.4 

0.01 

-.0000015 

.0005 

0.1 

2.18 

-.0004517 

0.15 

0.70 

-.0001432 

0.2 

0.30 

-.0000615 

0.3 

0.009 

-.0000181 

0.4 

0.04 

-.0000074 

.001 

0.1 

3.87 

-.0008032 

0.15 

1.31 

-.0002699 

0.2 

0.58 

-.0001189 

0.3 

0.18 

-.0000357 

0.4 

0.08 

-.0000147 

The  support  for  this  calculation  proceeds  as  follows:  Suppose  X  is  a 
random  variable  uniformly  distributed  on  the  interval  (-AH,AH).  Then  the 
mean  of  X  is  zero  and  the  variance  is  (AH)2/3.  If  there  are  n  layers  to  be 
processed  then  the  error  in  determining  the  horizontal  distance  is  the  sum  of 
N  independent  and  identical  such  X's.    It  will  have  mean  zero  and  standard 


deviation  AHVn/3;  zero  plus  or  minus  two  standard  deviations  could  be  a 
significant  amount,  especially  for  small  elevation  angles. 
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APPENDIX  E.  CONVERSION  OF  FOUR  TRANSIT  TIMES  TO  INPUTS  FOR 

THE  RAY  TRACING  ALGORITHM 

The  method  in  current  use  for  converting  the  four  transit  times  (ti, ...,  t4 ) 
into  an  azimuth  angle,  <J>C,  an  elevation  angle,  6C  and  tac  an  estimated  transit 
time  from  the  source  to  the  acoustic  center  is  outlined  and  critiqued  below. 

The  two  angles  <j>c  and  6C  are  generated  from  a  description  that  assumes  a 
constant  value  u  for  the  speed  of  sound  for  all  points  in  the  array.  Then  the 
concept  of  an  "apparent  position,"  (X,Y,Z)  relative  to  the  acoustic  center  for 
the  sources  in  a  constant  speed  median  is  utilized  for  purposes  of  estimating 
the  two  angles.  The  apparent  position  must  satisfy  the  system  of  equations 

(X  +  D/2)2+(Y-D/2)2  +  (Z-D/2)2  =  v2t2 
(X-D/2)2  +  (Y  +  D/2)2  +  (Z-D/2)2  =  v2t2 
(X-D/2)2  +  (Y-D/2)2  +  (Z  +  D/2)2  =  u2f2 
(X-D/2)2+(Y-D/2)2  +  (Z-D/2)2  =  u2f2 

This  is  a  system  of  four  equations  in  three  unknowns.  Unless  the  values  of  ti, 
...,  t4  are  singularly  coherent,  there  are  an  infinity  of  solutions  for  X,  Y,  Z. 

The  current  policy  is  to  obtain  a  unique  solution  by  subtracting  the  fourth 
equation  successively  from  each  of  the  other  three,  leaving  a  three  by  three 
system  remaining.   I.e., 

2XD=v2[t2-t2) 

2YD=  v2(t2-tl)  (E.2) 

2ZD=v2(tl-t24) 
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and  a  unique  solution.    Other  rationales  supporting  this  approach  can  be 
found  in  [5]. 

At  this  point  an  adjustment  is  made  under  the  name  of  a  direction  cosine 
correction  (DCC).  The  quantity  D/2  is  added  to  each  component  of  the 
apparent  position  (in  effect  making  the  position  relative  to  the  C-phone)  and 
the  length  of  the  new  vector  is  computed  and  divided  by  ut4.  Call  this 

DCC  =  ij{X  +  D I  if  +  (Y  +  D 1 2)2  +  (Z  +  D  /  2)2  /  vt4  (£.3) 

The  translated  values  are  normalized  by  DCC  prior  to  returning  the  origin  to 
the  acoustic  center: 

X±D/2 

c         DCC 
Yc  =  Y  +  D/2-D/2  (£.4) 

y±D/2      2 

c         DCC 

Next  come  the  tilt  corrections,  required  because  the  corrected  apparent 
position  is  in  the  coordinate  system  cs(ac),  see  Appendix  B.  This  is 
accomplished  by  applying  the  approximation  of  eq.  (B.4),  i.e.,  that  coefficient 
matrix  (instead  of  the  exact  one  (B.5))  to  the  vector  OQ,  Yc,  Zc):  The  Z  rotation 
can  also  be  applied  at  this  point,  i.e.,  multiplying  by  P3.  Let  us  call  the  result  of 
all  this  (X,Y,Z)  and  form  the  functions 

sin(0c)  =  Z  /  Vx2  +  Y2+Z2 

sin(0c)  =  Y/Vx2  +  Y2  (£.5) 

cos(0c)  =  X/Vx2  +  Y2 

and,  from  these,  the  angles  6C  and  <t>c  can  be  found. 
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Finally,  there  is  need  to  construct  a  transit  time  to  the  acoustic  center 
because  it  is  not  measured.  Calling  the  values  tac/  the  proportionality 
adjustment  with  t4  is  used. 


^  =  tjx2  +Y2  +  Z2  /  ^(X  +  D  /  2)2  +  (Y  +  D  I  if  +(Z  +  D  /2)2        (E.6) 
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APPENDIX  F.  ALTERNATIVE  INITIALIZATION  OF  RAY  TRACING 

Current  methodology  utilizes  the  four  measured  ray  transit  times  (ti,  t2, 
t3,  t4)  from  source  to  the  X,  Y  Z,  C  hydrophones  (respectively),  and  converts 
them  to  (|)C/  8C/  and  tac/  the  azimuth  and  elevation  angles  and  an  estimated 
transit  time  from  source  to  the  acoustic  center.  It  is  seen  in  Section  5  that 
there  can  be  considerable  error  in  this,  especially  for  <|)c.  Our  alternative  is  to 
shift  the  acoustic  center  to  the  C-phone  and  to  ignore  t3.  The  resulting  angles 
will  have  much  smaller  error,  and  the  transit  time  t4  is  used  directly. 

We  proceed  to  develop  the  ray  azimuth  (spherical  coordinate  longitude) 
and  elevation  (spherical  coordinate  latitude)  angles  from  times  ti,  t^  U,  in  the 
range  coordinate  system.  (See  Appendix  B  for  a  general  discussion  of 
coordinate  systems.)  The  conversion  of  those  times  is  influenced  by  the  array 
orientation. 

Let  Sj,cj  be  the  sine  and  cosine  of  the  ith  Euler  angle  in  the  orientation  of 
the  3-D  sensor  array  j;  i  =  1,  2,  3.  The  conversion  of  XTILT,  YTILT,  and  ZROT 
into  Euler  angles  is  explained  in  Appendix  B,  together  with  methodology  for 
locating  the  positions  of  the  phones  at  the  ends  of  the  array  arms. 

Having  moved  the  acoustic  center  to  the  C-phone,  we  are  using  the 
coordinate  system  cs(a).  Letting  Bj  represent  the  jth  column  of  the  matrix  B, 
eq.  (B.3),  the  locations  of  the  X,Y,Z  phones  are,  respectively, 

DBX       DB2       DB3 
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and,  of  course,  the  C-phone  is  at  the  origin.  These  locations  in  the  range 
coordinate  system  can  be  found  by  adding  the  location  of  the  C-phone  as 
indicated  in  Table  B-2. 

Typically  there  is  very  little  difference  in  the  depths  (3rd  component)  of 
the  X,Y,  and  C-phones.  Hence  there  is  small  variability  if  the  speed  of  sound 
at  these  three  depths  and  the  use  of  a  constant  common  value,  say  u,  is  more 
tenable  than  the  previous  use. 

Again  we  adopt  the  concept  of  an  apparent  position  (X,Y,Z)  of  the  sound 
source  and  the  first  two  direction  cosines  can  be  found  by  applying  the  law  of 
cosines  (see  the  figure)  for  j  =  1,  2. 


where   ti,  t2,  t4  are  the  signal  transit  times  to  the  X,  Y  and  C  phones, 
respectively. 

For  j  =  1,  the  cosine  of  the  marked  angle  is  X/ut4  and  by  the  law  of  cosines 
(v  ti)2  =  D2+(\)t4)2  -  2D  •  t4  cos^) 


It  follows  that 


D 


X=  j  +  (u  •  t4 -  u  •  ti)(u  •  t4  +  u  •  ti)/(2D) 
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and  a  like  equation  for  j  =  2,  leading  to 

Y  =  D/2  +  (\)t4  -  ut2)  (ut4  +  \)t2)  /  2D. 

The  third  direction  cosine  is  obtained  from 

cosOF3)  =  Vl-cos2(xP1)-cos2(Y2)  =  A/ 1  -  (X2  +  Y2)  /  u2t 


Next  we  rotate  these  values  into  alignment  with  the  test  range  coordinate 
system,  i.e.,  apply  the  matrix  B,  eq.  (B.3) 


(F.l) 


The  proposed  ray  elevation  (0)  and  azimuth  (<|)p)  angles  can  be  found  from 

(F.2) 


X" 

"X" 

yc 

=  B 

Y 

-ZC 

Z 

cos(ep)  =  zjyjx[+rf+zt 


and 


sm(<pp)  =  Yc/^X?+Y? 
cos(0F)  =  Xc/Vxc2  +  Yc: 


(F.3) 


Comment.  This  technique  was  designed  to  treat  the  speculated  shortcomings 
of  the  procedure  currently  in  use.  The  results  are  quite  successful,  especially 
in  reducing  azimuth  error.  It  has  the  curious  property  of  not  using 
information  provided  by  the  Z-phone.  It  seems  wasteful  not  to  use  this 
information.  Yet  any  system  that  does  use  it  must  be  required  to  perform  at 
least  as  well  as  this  one. 
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Appendix  G. 

SUBROUTINE  ERRCUR(H1  ,Z1  ,K, A2,XTILT,YTILT/ZROT,D/L,VEL/M/ 
*  THR,THC,THER,PHR,PHC,PHER,HC,HER,ZC,ZER,TEE,TIMER) 

C  08/10/90 

C       COMPUTES  THE  TRUE  ELEVATION  ANGLES  AND  THE  ESTIMATED  ELEVATION 

C  ANGLES  FOR  K  SOUND  SOURCE  DISTRIBUTED  EQUALLY  ON  A  CIRCLE  OF  RADIUS 

C  HI  AT  DEPTH  Zl  FOR  A  RECIEVING  ARRAY  AT  DEPTH  A2  BUT  OTHERWISE  AT  THE 

C  CENTER  OF  THE  CIRCLE.  THE  ACCOUSTIC  CENTER  IS  THE  GEOMETRIC  CENTER  OF 

C  THE  ARRAY  CUBE.  THE  CURRENTLY  USED  METHODOLOGY  IS  APPLIED. 

C 

p»**»»»»»***»»*»»»»»**»**»*»»»»»***»»*»»»*»»*»»**»»»»**»»»»»»»»»» 

C 

C  INPUTS: 

C  HI:  RADIUS  OF  CIRCLE  IN  FEET 

C  Zl:  DEPTH  OF  SOUND  SOURCE  IN  FEET 

C  K:  NUMBER  OF  SOURCES  ON  THE  CIRCLE 

C  A2:  DEPTH  OF  THE  C-HYDROPHONE 

C  XTILT,YTILT,ZROT:   ORIENTATION  INFORMATION  ABOUT  THE 

C  SENSING  ARRAY  (RADIANS) 

C  D  LENGTH  OF  ARRAY  EDGES. 

C  L:  DEPTH  OF  LAYER  BOUNDARIES. 

C  M:  NUMBER  OF  RECORDS  IN  THE  VELOCITY  DEPTH  PROFILE. 

C  VEL:  AVERAGE  SPEED  OF  SOUND  IN  THE  LAYERS. 

C 

C  OUTPUTS: 

C  THR:  ACTUAL  ELEVATION  ANGLE  (THETA) 

C  THC:  THETA  ESTIMATES  AT  ACCOUSTIC  CENTER,  CURRENT  METHOD 

C  THER:  THETA  ERROR  AT  THE  ACCOUSTIC  CENTER,  CURRENT  METHOD 

C  PHR:  ACTUAL  AZIMUTH  ANGLE  (PHI). 

C  PHC  PHI  ESTIMATES  AT  ACCOUSTIC  CENTER,  CURRENT  METHOD. 

C  PHER:  PHI  ERROR  AT  ACCOUSTIC  CENTER. 

C  HC  HORIZONTAL  ESTIMATE,  CURRENT  METHOD. 

C  HER:  HORIZONTAL  ERROR. 

C  ZC  VERTICAL  ESTIMATE,  CURRENT  METHOD. 

C  ZER:  VERTICAL  ERROR. 

C  T:  TRANSIT  TIME  TO  ACCOUSTIC  CENTER. 

C  TIMC:  ESTIMATE  OF  TIME  TO  ACCOUSTIC  CENTER,  CURRENT  METHOD. 

C  TIMER:   TRANSIT  TIME  ERROR. 

DIMENSION  B(5,3),DZ(300),THER(30),HD(5),HC(30),ZC(30) 
DIMENSION  L(300),LM(300),LL(300),PX(30),PY(30),T(5),TEE(30) 
DIMENSION  THC(30),THR(30),V0(300),V  1  (300),VEL(300),V V(300) 
DIMENSION  HX(5),HY(5),HER(30),ZER(30),TIMER(30) 
DIMENSION  X0(3),PHC(30),PHR(30),PHER(30),TIMC(30) 

RE AL»8  Al , A2, A2M, ANG,B,C1  ,C2,C3,D,DP,DW,DZ,THER,GG 
REAL»8H,H0,H1,HD,L,LM,LL,LPZ,PIE,P1,P2,PX,PY,S1,S2,S3 
RE AL»8  T,T0,THEC,TH0,TH1  ,THC,THR,HC,ZC,HER,ZER 
REAL»8  V,V0,Vl,VEL,W,X,X0,XTILT,Y,YTILT,Z,Z0,Zl,ZROT 
REAL*8CC1,CC2,VV1,VV0,HX,HY,SC,DCC,RAC,RC,TEE 
REAL*8CAZ,SAZ,PHC,PHR,PHER,LH,C,T3P,TIMC,TIMER 
REAL»8  SR,SRER,LB,SV,SVU2,SVU,SU2,SU4,G1,G 

PIE  =  3.14159265359D0 
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IEST  =  0 
M1  =  M 

C    DISTRIBUTE  THE  K  SOURCES  EQUALLY  AROUND  THE  CIRCLE  COUNTER 
C    CLOCKWISE  FROM  THE  EAST. 
DO  10 1  =  1,K 

ANG  =  2*PIE*(I-1)/K 
PX(I)  =  Hl*DCOS(ANC) 
PY(I)  =  Hl'DSIN(ANG) 
PHR(I)  =  ANG 

IF(ANG.GT.PIE)  PHR(I)  =  ANG  -  2TIE 
10     CONTINUE 

C    FORM  SINES  AND  COSINES  OF  ALL  THE  EULER  ANGLES:ROLL,PITCH,YAW 

52  =  DSIN(XTILT) 

C2  =  DSQRTd  -  S2**2) 
SI  =  DSIN(YTILT)/C2 
CI  =  DSQRTd  -  Sl**2) 

53  =  -DSIN(ZROT) 
C3  =  DCOS(ZROT) 

C    IN  THE  COORDINATE  SYSTEM  HAVING  CENTER  AT  THE  C-HYDROPHONE 
C    AND  POSITIVE-UPWARD,  THE  LOCATIONS  OF  THE  FOUR  HYDROPHONES 
C    (RELATTVE  TO  THE  ARM  LENGTH  D)  ARE  DEVELOPED  NEXT.  THIS  IS 
C    THE  TRANSPOSE  OF  THE  MATRIX  B  IN  APPENDIX  B. 

B(l,l)  =  C2*C3 

B(l,2)  =  C2*S3 

B(U)  =  S2 

B(2,l)  =  -S1*S2*C3  -  C1*S3 

B(2,2)  =  -S1*S2*S3  +  C1»C3 

B(2,3)  =  S1*C2 

B(3,l)  =  -C1*S2»C3  +  S1*S3 

B(3,2)  =  -C1*S2*S3  -  S1*C3 

B(3,3)  =  C1*C2 

C    LIKE  NOTATION  WILL  BE  USED  TO  LOCATE  THE  C-HYROPHONE  AND  THE 
C    ACCOUSTIC  CENTER  (AC). 
D012J  =  13 
B(4J)  =  0.0D0 

B(5,])  =  0.5*(B(1,J)  +  B(2,J)  +  B(3J)) 
12     CONTINUE 

Al  =  0.0D0 
P2  =  Z1 

C    LOCATE  THE  HYDROPHONE  HORIZONTAL  COMPONENTS  IN  THE  COORDINATE 
C    SYSTEM  CENTERED  AT  AC. 
DO  14  J  =  1,5 

HX(J)  =  D*(B(J,1)-B(5,U) 
HY(J)  =  D*(B(J,2)  -  B(5,2)) 
14     CONTINUE 

C  DETERMINE  THE  DEPTHS  OF  THE  FOUR  HYDROPHONES  AND  THE  AC. 
HD(1)  =  A2  +  D*(B(5,3)  -  B(l,3)) 
HD(2)  =  A2  +  D*(B(5,3)  -  B(2,3)) 
HD(3)  =  A2  +  D*(B(5,3)  -  B(3,3)) 
HD(4)  =  A2  +  D*(B(5,3)  -  B(4,3)) 
HD(5)  =  A2 
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C   FIND  THE  DEEPEST  HYDROPHONE 
A2M  =  0.D0 
D051J=1,4 

IF(HD(J).GT.A2M)  A2M  =  HD(J) 
51      CONTINUE 

C    FORM  THE  SET  OF  LAYER  MIDPOINTS. 

DO  105 1  =  1,M-1 

LM(I)  =  .5*(L(I)  +  LU+1)) 
105    CONTINUE 

LM(M)  =  LM(M-l)  +  L(M)  -  L(M-l) 

C      FORM  DEPTH  INCREMENTS,  AND  ALL  SOUND  VELOCITY  SLOPES  AND 
C     INTERCEPTS. 
DO110I=l,M-l 

DZ(I)=LM(I+1)-LM(I) 

V0(I)=(LM(I+1)»VEL(I)  -  LM(I)*VEL(I+1))/DZ(I) 
V1(I)=(VEL(I+1)-VEL(I))/DZ(I) 
110    CONTINUE 
C 

IF(A2M.LT.LM(M-2))  GOTO  126 

C  IF  A2M  IS  DEEPER  THAN  THE  LAST  LAYER  MIDPOINT,  THEN  WE  EXTRAPOLATE 
C  THE  SOUND  VELOCITY  PROFILE  BY  USING  A  QUADRATIC  FUNCTION  OVER  THE 
C   DEEPEST  100  FEET. 

C    FIRST  COUNT  THE  NUMBER  OF  LAYERS  (OF  THICKNESS  DZ(M-2))  TO 
C    BE  ADJOINED.  ALSO  MUST  EXTEND  THE  L  ARRAY. 

K0  =  2  +  MAX(0,NINT((A2M-LM(M-1))/DZ(M-1))) 

MC  =  21 

C    FIND  THE  AVERAGE  DEPTH  FOR  THE  LAST  100  FEET. 

LB  =  0.0D0 

DO200J  =  M+l-MC,M 
200  LB  =  LB  +  LM(J) 

LB  =  LB/MC 

C    FORM  SUM  OF  POWERS  AND  PRODUCTS. 
SV  =  0.0D0 
SVU2  =  0.0D0 
SVU  =  0.0D0 
SU2  =  0.0DO 
SU4  =  0.0D0 
Gl  =  0.0D0 
DO210J  =  M+l-MC,M 

U  =  LM(J)  -  LB 

SV  =  SV  +  VEL(J) 

SVU  =  SVU  +  U»VEL(J) 

SVU2  =  SVU2  +  U**2  *  VEL(J) 

SU2  =  SU2  +  U»*2 

SU4  =  SU4  +  U**4 

G1=G1  +  V1(J) 
210    CONTINUE 

G1=G1/MC 

G  =  SVU/SU2 

GG  =  (MOSVU2  -  SU2*SV)/(SU4*MC  -  SU2»»2) 
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IF(GG.LT.O.ODO)  GG  =  O.ODO 
IF(Vl(M-l).LT.O.ODO)  Vl(M-l)  =  Gl 

C    PERFORM  THE  EXTRAPOLATION. 
DO  125  I=M,M+K0 

V1(I)  =  Vl(I-l)  +  GG*DZ(M-1) 

LM(I+1)  =  LM(I)  +  DZ(M-l) 

VELd+1)  =  VEL(I)  +  DZ(M-1)*V1(I) 

V0(I)  =  (LM(I+1)*VEL(I)  -  LM(I)»VEL(I+1))/DZ(M-1) 

L(I+1)  =  L(I)  +  DZ(M-1) 

DZ(I)  =  DZ(M-1) 

125  CONTINUE 

C    UPDATE  M,  THE  NUMBER  OF  LAYERS 
M  =  M+KO 

126  CONTINUE 

C    THE  OUTER  LOOP  WILL  PERFORM  COMPUTATIONS  FOR  THE  K  SOUND 
C    SOURCES. 

C    ADJUST  DV  TABLE  TO  25  FT.  INCREMENTS. 
CALL  VELMOD(L,VEL,Ml,LL,VV,MM) 

C    RAYFITTING 

C    THE  INNER  LOOP  WILL  FIT  RAYS  TO  THE  FOUR  HYDROPHONES 
C    AND  THE  AC  IN  THE  ORDER  X,Y,Z,C  AND  AC. 
DO  50 1  =  1,K 

WRITE(V)'  OUTER  LOOP  I  =  ',1/  K  =  ',K 
D035J=L5 

PI  =  DSQRT((PX(I)-HX(J))**2  +  (PY(I)-HY(J))**2) 
Z0  =  HD(J) 

CALLRAYFIT1(A1,Z0,P1,P2,M,VEL,LM,DZ,V0,V1,T0,TH0, 
TH1JEST) 

C    COLLECT  THE  FIVE  TRANSIT  TIMES. 
T(J)  =  TO 

C    IN  THIS  PROGRAM  WE  KEEP  ONLY  THE  TRUE  ELEVATION  ANGLE  AT  AC. 

THR(I)  =  THO 
35  CONTINUE 

C   INNER  LOOP  COMPLETED. 

C    LOCATE  THE  WATER  LAYER,  N,  CONTAINING  THE  ARRAY.  USE  THIS 
C    TO  DEVELOP  THC,  THE  CURRENTLY  USED  ESTIMATE  OF  THO. 
N  =  MM 
DO  37  J  =  2,MM 

IF((LL(J-1).LE.A2).AND.(LL(J).GT.A2))  N  =  J-l 
37  CONTINUE 

C 

V  =  VV(N) 

C    USE  THE  FOUR  TRANSIT  TIMES  TO  PRODUCE  ESTIMATES  OF  THE 

C    ENTRANCE  ANGLE.  CALCULATE  THE  PRE-TILT  CORRECTED  APPARENT 

C    POSITION  AND  INCLUDE  THE  DIRECTION  COSINE  CORRECTION. 

DO40J  =  U 

X0(J)=((V*T(4)-V*T(J))*(V*T(4)+V*T(J)))/(2*D) 

X0(J)  =  0.5*D  +  X0(J) 
40  CONTINUE 


SC  =  V*T(4) 

DCC  =  (DSQRT(X0(1)*»2  +  X0(2)»»2  +  X0(3)»*2))/SC 

D041J=U 

X(XJ)  =  (XO(J)/DCC)  -  0.5*D 
41  CONTINUE 

C    NEXT  MAKE  TILT  CORRECTIONS 
X  =  X0(1)  -  X(X3)»DSIN(XTILT) 
Y  =  X0(2)  -  X(K3)»DSIN(YTILT) 
Z  =  X0(3)  +  X0(1)*DSIN(XTILT)  +  X0(2)»DSIN(YTILT) 

RAC  =  DSQRT(X*»2  +  Y»*2  +  Z»»2) 

RC  =  DSQRT((X  +  D/2)*»2  +  (Y  +  D/2)**2  +  (Z  +  D/2)*»2) 

TO  =  T(4)*RAC/RC 

TIMC(I)  =  TO 

C    PERFORM  Z  ROTATION  IN  THE  (X,Y)  PLANE. 
X(KD  =  C3»X-S3*Y 
X0(2)  =  S3»X  +  C3*Y 

C    COMPUTE  THEC:  THE  ESTIMATE  OF  THETA  CURRENTLY  IN  USE. 
THEC  =  DASIN(Z/DSQRT(X**2  +  Y**2  +  Z**2)) 

C    COMPUTE  SINES  AND  COSINES  OF  AZIMUTH. 
SAZ  =  X0(2)/DSQRT(X**2  +  Y**2) 
CAZ  =  X0(1)/DSQRT(X**2  +  Y»»2) 
PHC(I)  =  DATAN2(SAZ,CAZ) 
PHER(I)  =  PHC(I)  -  PHR(I) 

IF(ABS(PHER(I)).GT.PIE)  PHER(I)  =  PHC(I)  +  PHR(I) 
IFG  =  0 

C    RAYTRACE  BY  THE  ISOSPEED  METHOD. 

CALLISOSPEED(A1^0,TO,THEC,LL,VV/MM/H^/TH1) 

HC(I)  =  H 
ZC(I)  =  Z 

C    NOW  FINISH  THE  OUTPUT. 
THC(I)  =  THEC 

C    AND  THE  ERRORS 

TEE(I)  =  T(5) 

THER(I)  =  THC(I)  -  THR(I) 

TIMER(I)  =  TIMC(I)  -  T(5) 

HER(I)  =  HC(I)  -  HI 

ZER(I)  =  ZC(I)-Z1 

SR  =  DSQRT(H1**2  +  (A2  -  Zl)»*2) 

SRER  =  DSQRT(HC(I)**2  +  (A2  -  Zl)**2)  -  SR 
50     CONTINUE 
C   OUTER  LOOP  COMPLETED! 

RETURN 

100    FORMAT(3(5X,E13.6)) 
120    FORMAT(3(5X,F15.12)) 
130    FORMAT(3(F12.8,2X)) 
END 
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C    LIBRARY  FILE:  LIB12.FOR  2/22/90 

£*      *  +  **  +  **  +  ***  +  ***■***»**+*+***  +  **+**  +  *  +  **»  +  *+****  +  ***■*•*  ++*++++  +  *  +  * 

SUBROUTINE  ISOGRAD1(A1/A2,TO/THO,N,LM,VEL/VO,V1,DZ,H/Z,TH1) 

C  09/25/89 

C  TO  :  TRANSIT  TIME  (SEC). 

C  THO :  ELEVATION  ANGLE  AT  SENSOR  (RAD). 

C  Al  :  HORIZONTAL  COORDINATE  OF  SENSOR. 

C  A2  :  VERTICAL  COORDINATE  OF  SENSOR,  POSITIVE  DOWN. 

C  V0,V1  :  ARRAYS  CONTAINING  SOUND  VELOCITY  PARAMETERS. 

C  LM  :  ARRAY  CONTAINING  LAYER  MIDPOINTS. 

C  N  :  INDEX  OF  DEEPEST  LAYER  USED. 

C 

DIMENSION  LM(300),V0(3OO),V  1  (300),DZ(300),VEL(300) 

REAL»8T0,H,H0,Z,A1,A2,TH0,TH1,VEL 

REAL*8  LM,V0,V1,DZ,Q1,Q2 

REAL*8  VA2,R,VP2,TH,RV,DW,DT,X,T 

INTEGER  N,IS 

I  =  N 

T  =  0.0D0 

TH=TH0 

H0=A1 

VA2=V0(I)+V1(I)*A2 

RV=DCOS(TH)/VA2 

Z  =  A2 

DZ(N)  =  Z  -  LM(N) 

50  IF(VKI).EQ.O.O)  THEN 

DW  =  DZ(I)/DSIN(TH) 
DT  =  DW/VOU) 
H  =  HO  +  DWDCOS(TH) 
TH1=TH 
ELSE 

Q2=-V(XI)/V1(I) 
IF  (Q2)  51,52,53 

51  IS  =  -1 
GOTO  54 

52  IS  =  0 
GOTO  54 

53  IS=1 

54  CONTINUE 

Q1=H0  +  (Q2-Z)*DTAN(TH) 

R=DSQRT((Q2-Z)**2  +  (Ql-H0)**2) 

THl=DACOS(RV*VEL(D) 

DT=DLOG((DCOS(TH)/(l+DSIN(TH)))*((l+DSIN(THl))/ 
DCOS(THl)))/Vl(I) 

H=Q1  -  IS*R*DSIN(TH1) 
ENDIF 
T=T+DT 

IF(T.GE.T0)GOTO60 
Z=LM(I) 
H0  =  H 
TH=TH1 
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1=1-1 

GOTO  50 

60     DT=T0+DT-T 

IF(V1(I).EQ.0.0)  THEN 

DW  =  V0(I)*DT 

DZ(I)  =  DWDSIN(THl) 

H  =  HO  +  DWDCOS(THl) 

Z  =  Z-DZ(I) 
ELSE 

X=(EXP(DT»Vl(I)))*(l+DSIN(TH))/DCOS(TH) 

THl=DACOS((2»X)/(l+X*»2)) 

H  =  Ql  -  IS»R*DSIN(TH1) 

Z  =  Q2  -  IS*R*DCOS(THl) 
ENDIF 

C    RESTORE  THE  END  LAYERS. 

DZ(I)  =  LM(I+1)-LM(I) 
DZ(N)  =  LM(N+1)  -  LM(N) 

RETURN 

END 

/-•******»*********»*»******•<■********************  +  **»********»**»** 


SUBROUTINE  ISOSPEED(A1,A2,TO,THO,L,VEL/M,H/Z,TH1) 

C  08/09/89 

C  This  is  a  2-D  ray  tracing  algorithm  that  mimics  the  one  in 

C  Proceedure  5181.  It  utilizes  the  assumption  that  the  speed 

C  of  sound  in  water  is  constant  for  the  entire  layer  encompassed 

C  by  the  layer  boundaries.  A  fixed  ray  invariant  is  used 

C     throughout  the  entire  migration. 

^************* **********************  ***************************** 

C  INPUTS: 

C  A1,A2  -  POSITION  OF  SENSOR  (A2>0  DOWN) 

C  TO  -  TRANSIT  TIME 

C  THO  -  ELEVATION  ANGLE  (OF  THE  RAY  AT  THE  SENSOR, 

C  ALSO  CALLED  THE  ENTRANCE  ANGLE) 

C  L  -  ARRAY  CONTAINING  LAYER  BOUNDARIES 

C  VEL  -  ARRAY  CONTAINING  SOUND  VELOCITY  AT  THE 

C  THE  LAYER  MIDPOINTS 

C  OUTPUTS: 

C  HZ  -  POSITION  OF  TARGET  (SOUND  SOURCE) 

C  TH1  -  ELEVATION  ANGLE  AT  TARGET 

DIMENSION  L(300),VEL(300) 

REAL»8  Al,A2,C,SXDT,DWJH,DZ,RV/rHl,Z,H,L,VEL,TH0,T0 

Z  =  A2 

H  =  A1 

T  =  0.0 

C    CHOOSE  N  SUCH  THAT  L(N)  <=  A2  <  L(N+1).  IF  A2  IS  DEEPER 

C    THAN  LOWEST  LAYER  THEN  N  =  M,  THE  INDEX  OF  THE  DEEPEST  LAYER 

C    BOUNDARY 

N  =  M 


83 


DO  5  I  =  2,M 

IF((L(I-1).LE.A2).AND.(L(I).GT.A2))  N  =  I  -  1 
5     CONTINUE 

RV  =  DCOS(TH0)/VEL(N) 

J  =  N 

TH  =  TH0 

S  =  DSIN(THO) 

C  =  DSQRTd  -  S**2) 
10    DZ  =  Z-L(J) 

C    COMPUTE  THE  INCREMENTAL  SLANT  RANGE 
DW  =  DZ/S 

C    COMPUTE  THE  INCREMENTAL  TRAVEL  TIME 
DT  =  DW/VELO) 

C    ACCUMULATE  TOTAL  TRAVEL  TIME  AND  TEST 
T  =  T  +  DT 
IF  (T.GE.TO)  GOTO  50 

C    UPDATE  THE  HORIZONTAL  AND  VERTICAL  ACCUMULATIONS 
H  =  H  +  DW*C 
Z  =  Z-DZ 

C    USE  SNELL'S  LAW  TO  UPDATE  THE  LAYER  ENTRANCE  ANGLE  AND 
C    THE  TRIG  FUNCTIONS. 

J  =  J-1 

C  =  RV  *  VEL(J) 

S  =  DSQRTd  -  C**2) 

GOTO  10 

50    T  =  T-DT 
DT  =  TO  -  T 
DW  =  VEL(J)*DT 
DZ  =  DW*S 
H  =  H  +  DW*C 
Z=Z-DZ 
TH1  =  DASIN(S) 

RETURN 
END 


SUBROUTINE  RAYFIT1(A1,A2,P1,P2/M,VEL,LM,DZ/V0,V1,T0/TH0, 
TH1,IEST) 

C  09/12/89 

C    NEW  SUBROUTINE  TO  REPLACE  TGEN,  RAYTRACING  ALGORITHM. 

C  INPUTS: 

C  A1,A2  -  POSITION  OF  SENSOR  (A2  >  0  DOWN) 

C  P1,P2  -  POSITION  OF  SOUND  SOURCE  (  P2  >  0  DOWN  ) 

C  LM  -  ARRAY  CONTAINING  LAYER  MIDPOINTS 

C  M  -  NUMBER  OF  LAYER  MIDPOINTS 

C  VEL  -  ARRAY  CONTAINING  SOUND  VELOCITY  AT  THE 

C  LAYER  MIDPOINTS. 
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C  VO         -  SPEED  INTERCEPT  VALUES 

C  VI  -  SPEED  SLOPE  VALUES 

C  DZ         -  DEPTH  INCREMENTS 

C  IEST  -  FLAG  FOR  INITIALIZING  THE  ANGLE 

C  OUTPUTS: 

C  TO  -  TRANSIT  TIME 

C  THO       -  ELEVATION  ANGLE  AT  THE  SENSOR 

C  TH1       -  ELEVATION  ANGLE  AT  THE  SOUND  SOURCE 

DOUBLE  PRECISION  VEL(300),DZ(300),LM(300),V0(300) 
DOUBLE  PRECISION  V1(300),ANG(300),G(300) 
REAL-8  A1,A2,P1/P2,T0/TH0/TH1,EP,S,C 
REAL»8  H,H0,DW,VA2,VP2,GG,R,Z/TH,RV/Q1,Q2 
INTEGER  M,IS 

EP  =  lD-6 

C    DETERMINE  LAYERS  INVOLVED  IN  RAY  FITTING 
N  =  M 
J  =  M 
DO30I=l,M-l 

IF  ((LM(I).LE.A2).AND.(LM(I+1).GT.A2))  N=I 
IF  ((LM(I).LE.P2).AND.(LM(I+1).GT.P2))  J=I 
30  CONTINUE 

C    MAKE  END  CORRECTIONS  FOR  THE  LAYERS 
DZ(N)  =  A2  -  LM(N) 
DZ(J)  =  LM(J+D-P2 

C       COMPUTE  SPEED  OF  SOUND  AT  A2  AND  P2 
VA2  =  V0(N)  +  V1(N)*A2 
VP2  =  V0(J)  +  V1(J)T2 
IF(IESTJME.O)  GOTO  50 

C    INITIALIZE  THE  ELEVATION  ANGLE  AT  THE  SENSOR,  THO,  BY 
C    FITTING  A  STRAIGHT  LINE  SPEED  PROFILE  BETWEEN  P2  AND  A2. 
IF(VEL(N).EQ.VEL(J))  THEN 

THO  =  DATAN((A2-P2)/(P1-AD) 
ELSE 

Q2  =  (VEL(N)»LM(J)  -  VEL(J)*LM(N))  /  (VEL(N)-VEL(J)) 
Ql  =  0.5*(P1+A1)+(0.5»(P2-A2)*(P2+A2-2»Q2))/(P1-A1) 
THO  =  DATAN((Q1-A1)/(Q2-A2)) 
ENDIF 

C    OUTER  LOOP:  SET  UP  RAY  FITTING  FOR  THO  =  ELEVATION  ANGLE 
50    S  =  DSIN(TH0) 

C  =  DSQRTU.O  -  S**2) 

I  =  N 

RV  =  C/VA2 

H0  =  A1 

Z  =  A2 

60    IF(V1(I).EQ.0.0)  THEN 
DW  =  DZ(I)/S 
H  =  HO  +  DW»C 
ELSE 

Q2  =  -V0(I)/V1(I) 
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IF  (Q2)  61,62,63 

61  IS  =  -1 
GOTO  64 

62  IS  =  0 
GOTO  64 

63  IS=1 

64  CONTINUE 

Ql  =  HO  +  (Q2-Z)*S/C 
R  =  DSQRT((Q2-Z)**2  +  (Ql-H0)**2) 
C  =  RV*VEL(I) 
S=  DSQRT(1.0-C**2) 
H  =  Ql  -  IS*R*S 
ENDIF 

IF(I.EQ.J)GOTO80 
H0  =  H 
Z  =  LM(I) 
1  =  1-1 
GOTO  60 

80    TH1  =  DACOS(RV»VP2) 

C    FRACTIONAL  LAYER  CORRECTION 

IF(V1(J).NE.0.0)  H  =  Ql  -  IS*R*DSIN(TH1) 
IF  (ABS(H-Pl).LT.EP)  GOTO  90 

C    RE-ESTIMATE  THO 

THO  =  DATAN(DTAN(TH0)*H/P1) 
GOTO  50 

C    PREPARE  FOR  COMPUTATION  OF  TRANSIT  TIME. 
C    COLLECT  EXIT  AND  ENTRANCE  ANGLES. 
90    ANG(J)  =  TH1 
ANG(N+1)  =  THO 
D0  95I  =  J+1,N 

ANG(I)  =  DACOS(RV*VEL(I)) 
95    CONTINUE 

C    COMPUTE  TRANSIT  TIME 
TO  =  0.0D0 
DO  100 1  =  J,N 

IF(VKI).EQ.O.O)  THEN 

TO  =  TO  +  DZ(I)/(V0(I)*DSIN(ANG(I))) 
ELSE 

T0=T0  +  DLOG((DCOS(ANG(I+l)),f(l+DSIN(ANG(I))))/ 
((l+DSIN(ANG(I+l)))*DCOS(ANG(I))))/VKI) 
ENDIF 
100    CONTINUE 

C    REMOVE  THE  END  CORRECTIONS. 
DZ(J)  =  LM(J+1)-LM(J) 
DZ(N)  =  LM(N+1)  -  LM(N) 
IEST  =  1 

RETURN 
END 

/"-*!(■*****   +   *»**»***»**«.**»**»***»   +   **********»*******»*»*»**»**  +   +  **» 
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SUBROUTINE  VELMOD(L,VEL,M,LL,VV,MM) 

C  02/22/90 

C  THIS  PROGRAM  TAKES  THE  VELOCITY  SOUND  PROFILE  GIVEN  IN  FTVE  (5) 

C  FOOT  INCREMENTS  AND  CONVERTS  IT  INTO  TWENTYFIVE  (25)  FOOT  INCREMENT 

C  PROFILE. 

C 

C  INPUT: 

C        L:  DEPTH  IN  5  FT  INCREMENTS 

C        VEL:      SOUND  VELOCITY  IN  5  FT.  INCREMENTS 

C        M:         NUMBER  OF  ELEMENTS  IN  DEPTH  ARRAY 

C  OUTPUT 

C        LL         DEPTH  IN  25  FT  INCREMENTS 

C        VV:        SOUND  VELOCITY  IN  25  FT  INCREMENTS 

C        MM:      NUMBER  OF  ELEMENTS  IN  DEPTH  ARRAY 

£»»»»»»»»»»»»»***»»»»»»*»»»#»»»*»»»»*»»»»»»»»»»»»*»**»»**»»»»»»»»*»»»» 

DIMENSION  L(300),LL(300),VEL(300),W(300) 
REAL»8  L,LL,VEL,VV,VS 

MM  =  INT(M/5) 
MREM  =  M  -  5»MM 
VS  =  0.0D0 
DO10J  =  l,MM 

LL(J)  =  L(5*J-4) 

W(J)  =  0.2*(VEL(5*J-4)  +  VEL(5*J-3)  +  VEL(5*J-2)  + 
VEL(5*J-1)  +  VEL(5*])) 
10     CONTINUE 

DO20J  =  l,MREM 

VS  =  VS  +  VEL(M+1-J) 
20     CONTINUE 

LL(MM+1)  =  L(5*MM  +  1) 
VV(MM+1)  =  VS/MREM 
RETURN 
END 


SUBROUTINE  ERRPROP(Hl,Zl ,K/A2,XTILT/YTILT/ZROT,D,L,VEL/M/ 
•    THR/THCER/THl  ER,PHR,PHER,PH1  ER,H2ER,Z2ER,T4) 

C  08/22/90 

C  POSITION  ERROR  ANALYSIS  WHEN  THE  ORIGIN  IS  OVER  THE  C-HYDROPHONE 

C  AND  THE  PROPOSED  SYSTEM  IS  APPLIED.  COMPUTES  THE  TRUE  ELEVATION 

C  ANGLES  AND  ESTIMATED  ELEVATION  ANGLES  FOR  K  SOUND  SOURCES 

C  DISTRIBUTED  EQUALLY  ON  A  CIRCLE  OF  RADIUS  HI  AT  DEPTH  Zl  FOR  A 

C  RECIEVING  ARRAY  AT  DEPTH  A2  BUT  OTHERWISE  AT  THE  CENTER  OF  THE 

C  CIRCLE.  A2  IS  DEPTH  OF  THE  GEOMETRIC  CENTER  OF  THE  ARRAY  CUBE.  THE 

C  ACCOUSTIC  CENTER  IS  THE  C-HYDROPHONE.  BECAUSE  OF  DIRECT 

C  MEASUREMENT  THERE  IS  NO  TRANSIT  TIME  ERROR. 

C 

c 

C  INPUTS: 

C  HI:        RADIUS  OF  CIRCLE  IN  FEET 
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C  Zl:         DEPTH  OF  SOUND  SOURCE  IN  FEET 

C  K:  NUMBER  OF  SOURCES  ON  THE  CIRCLE 

C  A2:         DEPTH  OF  THE  CENTER  OF  THE  ARRAY 

C  XTILT,YTILT,ZROT:   ORIENTATION  INFORMATION  ABOUT  THE 

C  SENSING  ARRAY  (RADIANS) 

C  D  LENGTH  OF  ARRAY  EDGES. 

C  L-  DEPTH  OF  LAYER  BOUNDARIES. 

C  M:  NUMBER  OF  RECORDS  IN  THE  VELOCITY  DEPTH  PROFILE. 

C  VEL:      AVERAGE  SPEED  OF  SOUND  IN  THE  LAYERS. 

C 

C  OUTPUTS: 

C  THR:         ACTUAL  ELEVATION  ANGLE  (THETA) 

C  THONE:  THETA  ESTIMATE,  PROPOSED  METHOD 

C  THER:  THETA  ERROR,  PROPOSED  METHOD 

C  PHR:         ACTUAL  AZIMUTH  ANGLE  (PHI). 

C  PHI:  PHI  ESTIMATES,  PROPOSED  METHOD. 

C  PH1ER:  PHI  ERROR,  PROPOSED  METHOD. 

C  HC:  HORIZONTAL  ESTIMATE,  CURRENT  METHOD. 

C  HER:         HORIZONTAL  ERROR. 

C  ZC  VERTICAL  ESTIMATE,  CURRENT  METHOD. 

C  ZER:         VERTICAL  ERROR. 

C  T:  TRANSIT  TIME  TO  THE  C-HYDROPHONE. 

C 

DIMENSION  B(5,3),DZ(300),HD(5),HX(5),HY(5),ITH2(20),L(300) 
DIMENSION  LM(300),PH1(20),PH1ER(20),PHC(20),PHER(20) 
DIMENSION  PHR(20),PX(20),PY(20),T(4),TH1ER(20),TH2(20) 
DIMENSION  TH2ER(20),THC(20),THCER(20),THONE(20),THR(20) 
DIMENSION  V0(300),V1(300),VEL(300),X0(30),PZ(2,2),PC(2) 
DIMENSION  PP(2,3),E(3),TH3(20),TH3ER(20),TT3(3) 
DIMENSION  IFL(2),SSC(2),SP(2),H2(20),Z2(20),H2ER(20) 
DIMENSION  Z2ER(20),T4(30) 

REAL*8B,DZ,HD,HX,HY,L/LM,PH1,PH1ER,PHC,PHER,PHR,PX 
REAL»8  PY,T,TH1  ER,TH2,TH2ER,THC,THCER,THONE,THR,VO,V1 
REAL*8  VEL,X0,PONE,T4 

REAL*8A1,A2,A2M,ALPH1,ANG,C,C1,C2,C3,CAZ,CT/CX,CX0,CY,CY0 
REAL*8CZ,CZ0/D,DR/DR1,DP1,DP2,DTDS,EPS,F,FT,GG,H1,HH1,P1 
RE AL*8  P2,PIE,Q1  ,Q2,QI P,R0,R1  ,S,S0,S1  ,S2,S3,S AZ,ST,SS,T0,T3 
REAL*8  T3P,TH0,TH1,THEC,V,VV0,VV1,X,XTILT,Y,Y0,Y1,YTILT 
REAL»8  Z,ZO,Z1,ZROT,ZZI,SC,DT 

REAL*8PP,E,DEL,A,TH3,TH3ER,TT3,SSC,SP,PZ,PC 
REAL-8  A11,A12,A21,A22,B1,B2,H2,H2ER,Z2,Z2ER,THE 
REAL*8  LB,SV,SVU,SU2,SU4,SVU2,G,U 
INTEGER  ITH2,IFL 

PIE  =  3.141 59265359D0 
EPS=lD-6 
IEST  =  0 
M1  =  M 

C    DISTRIBUTE  THE  K  SOURCES  EQUALLY  AROUND  THE  CIRCLE  COUNTER- 
C    CLOCKWISE  FROM  THE  EAST. 
DO10I  =  l,K 

ANG  =  2*PIE*(I-1)/K 
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PX(I)  =  Hl*DCOS(ANG) 
PY(I)  =  H1*DSIN(ANG) 
PHR(I)  =  ANG 

IF(ANG.GT.PIE)  PHR(I)  =  ANG  -  2TIE 
10     CONTINUE 

C    FORM  SINES  AND  COSINES  OF  ALL  THE  EULER  ANGLES:ROLL,PITCH,YAW 

52  =  DSIN(XTILT) 

C2  =  DSQRTd  -  S2»»2) 
SI  =  DSIN(YTILT)/C2 
CI  =  DSQRTd  -  Sl*»2) 

53  =  -DSIN(ZROT) 
C3  =  DCOS(ZROT) 

C  IN  THE  COORDINATE  SYSTEM  HAVING  CENTER  AT  THE  C-HYDROPHONE 
C  AND  POSITIVE-UPWARD,  THE  LOCATIONS  OF  THE  FOUR  HYDROPHONES 
C    (RELATIVE  TO  THE  ARM  LENGTH  D)  ARE  DEVELOPED  NEXT. 

B(l,l)  =  C2*C3 

B(l,2)  =  C2»S3 

B(U)  =  S2 

B(2,l)  =  -S1*S2*C3  -  C1*S3 

B(2,2)  =  -S1*S2*S3  +  C1*C3 

B(2,3)  =  S1*C2 

B(3,l)  =  -C1*S2*C3  +  S1»S3 

B(3,2)  =  -C1»S2*S3  -  S1*C3 

B(3,3)  =  C1*C2 

C  LIKE  NOTATION  WILL  BE  USED  TO  LOCATE  THE  C-HYROPHONE  AND  THE 
C    ARRAY  CENTER. 

D012J  =  U 
B(4J)  =  0.0D0 

B(5J)  =  0.5*(B(J,1)  +  B(J,2)  +  B(J3» 
12     CONTINUE 

Al  =  0.0D0 
P2  =  Z1 

C    LOCATE  THE  HYDROPHONE  HORIZONTAL  COMPONENTS  IN  THE  COORDINATE 
C    SYSTEM  CENTERED  AT  C-HYDROPHONE. 
DO  14  J  =  1,5 

HX(J)  =  D-B(J,1) 
HY(J)  =  D*B(J,2) 
14     CONTINUE 

C  DETERMINE  THE  DEPTHS  OF  THE  FOUR  HYDROPHONES  AND  THE  ARRAY  CENTER. 
HD(1)  =  A2  +  D*(B(5,3)  -  B(l,3)) 
HD(2)  =  A2  +  D*(B(5,3)  -  B(2,3)) 
HD(3)  =  A2  +  D*(B(5,3)  -  B(3,3)) 
HD(4)  =  A2  +  D*(B(5,3)  -  B(4,3)) 
HD(5)  =  A2 

C   FIND  THE  DEEPEST  HYDROPHONE 
A2M  =  0.D0 
D051J=1,4 

IF(HD(J).GT.A2M)  A2M  =  HD(J) 
51     CONTINUE 

C    FORM  THE  SET  OF  LAYER  MIDPOINTS. 

&9 


DO  105 1  =  1,M-1 

LM(I)  =  .5*(L(I)  +  L(I+1)) 
105    CONTINUE 

C      FORM  DEPTH  INCREMENTS,  AND  ALL  SOUND  VELOCITY  SLOPES  AND 
C      INTERCEPTS. 
DO110I=l,M-2 

DZ(I)=LM(I+1)-LM(I) 

V0(I)=(LM(I+1)»VEL(I)  -  LM(I)*VEL(I+1))/DZ(I) 
V1(I)=(VEL(I+1)-VEL(I))/DZ(I) 
110    CONTINUE 

LM(M)  =  LM(M-l)  +  DZ(M-2) 
C 

IF(A2M.LT.LM(M-D)  GOTO  126 

C    IF  A2M  IS  DEEPER  THAN  THE  LAST  LAYER  MIDPOINT,  THEN  WE  EXTRAPOLATE 
C    THE  SOUND  VELOCITY  PROFILE  BY  USING  A  QUADRATIC  FUNCTION  OVER 
C   THE  DEEPEST  100  FEET. 

C    FIRST  COUNT  THE  NUMBER  OF  LAYERS  (OF  THICKNESS  DZ(M-2))  TO 
C    BE  ADJOINED.  ALSO  MUST  EXTEND  THE  L  ARRAY. 

K0  =  2  +  MAX(0,NINT((A2-LM(M-l))/DZ(M-2))) 

C    FIND  AVERAGE  DEPTH  OF  LAST  100  FEET. 
LB  =  0.0D0 
DO  43 1  =  M-21,M-1 
LB  =  LB  +  LM(I) 
43     CONTINUE 
LB  =  LB/21 

C    FORM  SUMS  OF  POWERS  AND  PRODUCTS. 
SV  =  0.0D0 
SVU2  =  0.0D0 
SVU  =  0.0D0 
SU2  =  0.0D0 
SU4  =  0.0D0 
D045I  =  M-21,M-1 

U  =  LM(I)  -  LB 

SV  =  SV  +  VEL(I) 

SVU  =  SVU  +  U*VEL(I) 

SVU2  =  SVU2  +  U**2  *  VEL(I) 

SU2  =  SU2  +  U**2 

SU4  =  SU4  +  U**4 
45     CONTINUE 

G  =  SVU/SU2 

GG  =  (21*SVU2  -  SU2*SV)/(SU4  -  SU2»*2) 

V1(M-1)  =  G 

C    PERFORM  THE  EXTRAPOLATION. 
DO125I=M,M+K0 

Vl(I-l)  =  VKI-2)  +  GG*DZ(M-1) 

LM(I)  =  LM(I-l)  +  DZ(M-2) 

VEL(I)  =  VEL(I-l)  +  DZ(M-2)*V1(I-1) 

V0(I-1)  =  (LM(I)*VEL(I-1)  -  LM(I-l)*VEL(I))/DZ(M-2) 

L(I+l)  =  L(I)  +  DZ(M-2) 

DZ(I-l)  =  DZ(M-2) 
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125  CONTINUE 

C    UPDATE  M,  THE  NUMBER  OF  LAYERS 
M  =  M+KO 

126  CONTINUE 

C    LOCATE  THE  WATER  LAYER,  N,  CONTAINING  THE  ARRAY. 
N  =  M 
D0  37J  =  2,M 

IF((LM0-DLEHD(4)).AND.(LM0).GT.HD(4)))  N  =  J-l 
37     CONTINUE 
C 

V  =  V0(N)  +  V1(N)*HD(4) 

C    THE  OUTER  LOOP  WILL  PERFORM  COMPUTATIONS  FOR  THE  K  SOUND 
C    SOURCES. 

C    RAYFITTING 

C    THE  INNER  LOOP  WILL  FIT  RAYS  TO  THE  FOUR  HYDROPHONES 
C    IN  THE  ORDER  X,Y,Z,  AND  C. 
DO50I  =  l,K 

WRITE(V)'  OUTER  LOOP  I  =  ',1/  K  =  ',K 
IW1  =0 
LX)35J  =  1,4 

PI  =  DSQRT((PX(I)-HX(J))*»2  +  (PY(I)-HY(J))»*2) 
Z0  =  HD(J) 

CALLRAYnTl(Al,Z0,Pl,P2,M,VEL,LM,DZ,V0,Vl,T0/TH0, 
»  TH1,IEST) 

C    COLLECT  THE  FOUR  TRANSIT  TIMES. 
T(J)  =  TO 

C    IN  THIS  PROGRAM  WE  KEEP  ONLY  THE  TRUE  ELEVATION  ANGLE  AT  THE 
C  C-HYDROPHONE. 

THR(I)  =  THO 

T4(I)  =  T(4) 
35  CONTINUE 

C    CALCULATE  THE  PRE-TILT  CORRECTED  APPARENT  POSITION 
DO40J  =  U 

X0(J)=(D»»2  +  (V*T(4)-V»T(J))»(V»T(4)+V»T(J)))/(2»D) 
40  CONTINUE 

C   COMPUTE  DIRECTION  COSINES  . 
CX0  =  X0(1)/(V*T(4)) 
CYO  =  X0(2)/(V*T(4)) 
CZO  =  DSQRTd  -  CX0**2  -  CY0*»2) 

C    PERFORM  EXACT  TILT  CORRECTIONS  AND  THE  ROTATIONAL  ALIGNMENT. 
42  CX  =  B(1,1)»CX0  +  B(2,1)*CY0  +  B(3,1)»CZ0 

CY  =  B(1,2)*CX0  +  B(2,2)»CY0  +  B(3,2)»CZ0 
CZ  =  B(1,3)*CX0  +  B(2,3)»CY0  +  B(33)*CZ0 
IF(IWl.EQ.l)  THEN 

TH2(I)  =  0.5*PIE  -  DACOS(CZ) 

GOTO  49 
ENDIF 

SAZ  =  CY/DSQRT(CX*»2  +  CY**2) 
CAZ  =  CX/DSQRT(CX*»2  +  CY"2) 
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PHI  (I)  =  DATAN2(SAZ,CAZ) 
PHIER(I)  =  PHI  (I)  -  PHR(I) 

IF(ABS(PH1ER(I)).GT.PIE)  PHIER(I)  =  PHl(I)  +  PHR(I) 
THONE(I)  =  0.5TIE  -  DACOS(CZ) 
THIER(I)  =  THONE(I)  -  THR(I) 

49  CALL  ISOGRAD1(A1/ZO,TO,THONE(I),N/LM/VEL,VO,V1,DZ,H2(I), 

Z2(I),THE) 
H2ER(I)  =  H2(I)  -  HI 
Z2ER(I)  =  Z2(I)-Z1 

50  CONTINUE 

C   OUTER  LOOP  COMPLETED! 

RETURN 

100    FORMAT(3(5X,E13.6)) 
120    FORMAT(3(5X,F15.12)) 

130    FORM  AT(4(F1 0.8,2X),F13.8,1  X,F1 2.8,2X,2(F1 2.8,2X),F1 2.8) 
140    FORMAT(5X,'The  transit  time  to  the  z-phone  is  not  bracketed') 
END 
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